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ABSTRACT 


In recent military conflicts, ballistic missiles have been used to achieve military 
and psychological objectives. With the proliferation of weapons of mass destruction 
(WMD), the growing threat of ballistic missiles being used as a delivery platform for 
WMD by rogue nations or militant groups becomes a concern for many countries. 
Defense against such threats becomes increasingly important. 

There are different guidance laws for the missile interception of aerial targets. 
These include pursuit, proportional navigation (PN) guidance as well as its variants. A 
new guidance algorithm was developed by John A. Lukacs IV and Prof Yakimenko in 
2006 to intercept a ballistic missile during the boost phase by a missile interceptor. This 
TS guidance algorithm uses the direct method of calculus of variations that maximizes 
the kinetic energy transfer from a surface-launched missile to a ballistic missile target. 

A trade-off study was conducted by applying this guidance law in simulated 
ballistic missile interception. This study examines the interactions and trade-offs between 
the various critical parameters in the intercept solution, like the endgame intercept 
geometry, time-to-intercept and intercept altitude. It provides insights into the feasibility 
and limitations of the TS guidance algorithm. A literature review of the drag model used 
in the algorithm and comparison of the new guidance with the compensated PN guidance 
was also conducted. A new induced drag model was developed for future studies. 

The results verified that the trajectory-shaping guidance is feasible for the 
interception of ballistic missiles in the boost phase for a wide range of interceptor launch 
locations with respect to a ballistic missile detection point. A better understanding of the 
trade-offs between the key parameters allows users to optimize the performance of this 
guidance. 
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I. INTRODUCTION 


A. BACKGROUND 

In recent conflicts, notably the first Gulf War in 1992 and the Israeli-Lebanon 
conflict in 2007, the use of tactical ballistic missiles to achieve military and psychological 
objectives had alerted political and military leaders around the world of the immense 
threat posed by such weapons. With the proliferation of weapons of mass destruction 
(WMD), there has been a greater concern that tactical ballistic missiles will be used by 
rogue nations or militant groups, as a delivery platform for such weapons [1], The 
outcome can be devastating, with the loss of many innocent lives and large-scale 
destruction. Many countries have invested substantially in developing ballistic missile 
defense systems as a strategy against such threats. In particular, the United States 
Department of Defense (DoD) was directed to develop a conceptual framework in 
2002—the Ballistic Missile Defense System (BMDS), in which active defense of 
intercepting incoming ballistic missiles in all phases of its trajectory, is a major strategy 
[2], An integrated ballistic missile defense system was subsequently developed based on 
a layered, defense-in-depth strategy. In order to intercept a ballistic missile, flying at 
supersonic speed and high altitude, the interceptor must have the required response time, 
speed and accuracy. The interception can take place in any of the three distinct phases— 
boost, midcourse or terminal, as illustrated in Figure 1. 
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Figure 1. Phases of Ballistic Missile Trajectory [3] 

During the boost phase, the ballistic missile rocket motors are firing in order to 
accelerate the missile to a high trajectory altitude. The advantages of boost phase 
interception are that the hot and bright rocket exhaust plume makes detection and 
targeting easier, and decoys cannot be used during this phase. The disadvantages lie in 
the geographical siting of the interceptor system (which has to be close to hostile 
territory) and the short time to intercept, typically about 180 seconds. Figure 2 shows 
images of ballistic missile interception by ground-launched intercept missile. In some 
literatures, defense analysts and scientists believed that interception of ballistic missile by 
surface-launched interceptor using established guidance laws is not feasible. However, 
the trajectory-shaping guidance developed by Lukacs and Prof Yakimenko in 2006 [4] 
shows that boost phase interception is possible, with promising results on the 
effectiveness of such interception, subject to limitations. 
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Figure 2. Ballistic Missile Interception by Ground-launched Missile [5] 

In the mid-course phase, the ballistic missile is in space after the rocket bums out. 
The coast period through space before re-entering the atmosphere can take several 
minutes, up to 20 minutes for a long-range Inter-Continental Ballistic Missile (ICBM). 
The advantage of intercepting a ballistic missile in this phase is that there will be 
sufficient decision and intercept time as well as a greater flexibility in the geographical 
defensive position of the interceptor system. However, such interception would require a 
larger and, hence, heavier interceptor missile, complemented by sophisticated radar and 
other sensors to handle potential space-based decoys. 

In order to intercept a ballistic missile in the terminal phase, the interceptor 
missile would have to do so after the ballistic missile re-enters the atmosphere. The 
advantages lie in the requirement for smaller and lighter interceptor missile, less 
sophisticated radar and lower possibility of decoy as they are not likely to work in this 
phase. The disadvantages are the very short reaction time (in the region of 30 seconds or 
less), less defended geographical coverage and possible effect of hazardous materials 
over the target area in the case of detonation of chemical, biological or nuclear 
warhead(s) mid-air. 

This paper focuses on the interception of ballistic missiles in the boost phase. 
Intercepting a ballistic missile in its boost phase is the ideal solution for missile defense, 
since the missile is most vulnerable during this phase of its flight. This is done usually 
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over the launch territory. There are, however, many challenges associated with this phase. 
The interceptor will have to contend with large acceleration rates, very short reaction 
time and requires reliable scanning and tracking [5,6]. There are several systems are 
under development for conducting boost phase interception, they include ground-based 
missiles, airborne lasers and space-based intercept missiles. 

This study only looks at surface-based missile interceptor. The missile guidance 
law is one of the most important factors that determines the feasibility and effectiveness 
of the intercept solution. If intercept is possible, then the next step is the interceptor’s 
capability to kill the target. With a very short engagement time, the interceptor missile 
needs a high speed and energy to reach the target. This imposes a limitation on the size of 
warhead. Earlier concepts and studies recognized that that a simple warhead effect is not 
sufficient to destroy an ICBM, hence the initiation of the development of hit-to-kill 
technologies [7,8], This concept relies on relatively smaller high speed missile, guided 
accurately to the target and transferring maximum kinetic energy to the ballistic missile 
for a kill. It suggests some form of geometry control during impact. 

The actual impact geometry is not an important parameter for most guidance laws 
implemented in intercept missiles which uses warhead effect (proximity fuze) as their 
main kill mechanism. However, for a hit-to-kill endgame condition which demands for 
maximum transfer of kinetic energy to the target missile, the intercept geometry 
commands a greater attention and has to be controlled as an input to the guidance law. 

The objective of this thesis is to conduct a trade-off study to evaluate the 
effectiveness of the hit-to-kill trajectory-shaping guidance through simulation, using the 
guidance algorithm developed by Lukacs and Prof Yakimenko in 2006. The code 
developed for the guidance law will “generate the interceptor’s entire flight path in order 
to minimize the distance traveled, minimize the time to intercept, and maximize kinetic 
energy transfer by controlling the interception geometry while providing near-optimal 
flight path to interception. This will be done by utilizing the direct method of calculus of 
variations combined with inverse dynamics theory to reverse engineer in real time, an 
optimal flight path using the missile’s onboard sensors and computers”. [9] 
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B. THESIS ORGANIZATION 


This thesis consists of six chapters. Chapter I provides a brief introduction and 
overview of ballistic missiles interception in the boost phase and outlines the objectives 
of this thesis. 

Chapter II provides a literature review of the models and guidance program, 
developed by Lukacs and Prof Yakimenko in 2006, based on the trajectory-shaping 
guidance law. The trade-off study is based on this guidance algorithm, with only slight 
modification in certain areas, for the simulation. For completeness, it is instructive to 
provide a brief overview of the guidance law, key characteristics and the MATLAB 
program developed as a prelude to the trade-off study and discussion. The detailed 
description of the theory behind the guidance algorithm and models used is appended in 
Appendix E. Much of the content in this chapter is extracted from the thesis written by 
Lukacs [10]. 

Chapter III discusses an alternative guidance law—the well-established 
Proportional Navigation (PN) guidance that is used in many advance missiles currently in 
service. It provides some background information on PN guidance and makes some 
comparison with the trajectory-shaping guidance. 

Chapter IV describes the drag model that is used in the trajectory-shaping 
guidance code and presents the need to incorporate an induced drag model as an 
enhancement to the existing model for better estimation of drag on both the ballistic and 
interceptor missile. 

Chapter V summarizes the results and discusses the feasibility of employing such 
guidance in a real-world scenario. 

The Appendices include data and plots of the simulation results obtained, report 
on the author’s participation in the AIAA Rocket Launch Competition, MATLAB code 
for the induced drag model, detailed description of the trajectory-shaping guidance and 
the complete MATLAB program codes. 
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II. SIMULATION AND MODELING 


A. DESCRIPTION OF TRAJECTORY SHAPING GUIDANCE 

In 2006, Lukacs and Prof Yakimenko [9] developed a simulation code in 
MATLAB as a demonstration of the feasibility of intercepting a ballistic missile during 
the boost phase by a surface-launched interceptor missile using the Trajectory Shaping 
guidance law (henceforth termed as TS guidance in the report). 

The TS guidance uses the principle of flight path optimization from the 
interceptor to the predicted target position. It relies on high-order polynomial as a 
reference function for the flight path and uses virtual as opposed to physical domain in 
optimization process. A preset thrust history is used in the computation of interceptor 
flight path. The flight path is derived by minimization a combined performance index 
including intercept geometry, time-to-intercept, penalty on altitude and dynamic 
constraints. The performance index, J is expressed in the following equation, 

J = tgo 2 + Wia(P - Pdesired) 2 + Pi + P2 

where t g0 is the time-to-intercept, p is the impact angle, Pi is the penalty on intercept 
altitude and P 2 is the penalty on dynamic constraints. The application of such guidance is 
particularly suitable for interception of targets with relatively fixed predicted trajectory, 
like that of a ballistic missile in the boost phase. 

In the simulation, the ballistic missile target was modeled after the Taepo-Dong 2 
(TPD-2) ballistic missile developed by the People’s Democratic Republic of Korea 
(PDRK, also known as North Korea) while the U.S.-made SM-6 Standard missile, an 
upgraded and extended range version of the SM-3 (shown in Figure 3) was used as the 
interceptor missile. A 3 degree-of-freedom (3DoF) mathematical model was previously 
developed and used to simulate the trajectory and flight characteristics of both the 
ballistic and interceptor missile based on available missile data from the open source. The 
intercept path is continuously calculated onboard the interceptor missile as a two-point 
boundary value problem, using Direct Methods of Calculus of Variation to calculate a 
near-optimal flight path and the control commands necessary to achieve it (refer to 
Appendix 3 for the detailed description of the TS guidance). 
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Figure 3. TPD-2 ICBM (left) and SM-3 Standard Missile (right) [11,12] 


The TPD-2 ballistic missile, developed by North Korea is believed to be an 
intercontinental ballistic missile (ICBM) capable of an effective range of 6000 to 6500 
km, with a 3-stage booster [13]. This put Alaska and all Asia Pacific countries up to 
northern Australia within reach of the TPD-2, from launch sites within North Korea [see 
Figure 4], In order to have a feasible solution of intercepting this missile during its boost 
phase, the interceptor must be fired from a ship or land-based launcher within range of 
the ballistic missile launch site. 



Figure 4. Reach of the TPD-2 ICBM Launched from North Korea 

The SM-6 (shown in Figure 5) is the U.S. Navy's latest Extended Range Anti-Air 
Warfare missile (ERAM) that employs active-homing terminal guidance—a key feature 
that allows for accuracy necessary to intercept an ICBM in the boost phase. The SM-6 
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missile is the upgraded and extended range version of the SM-3, currently under 
development by Raytheon, for the U.S. Navy as the next generation anti-ballistic missile 
defense system. It is designed to be capable of intercepting ballistic missiles in the boost 
phase if deployed within range of the ballistic missile launch site and complemented by 
suitable sensor systems to detect the launch. Figure 6 shows the comparison of the size of 
the TPD-2 and SM-6. 

M—■.^ 

Figure 5. The SM-6 Standard Missile [14] 



Figure 6. Comparison of the size of TPD-2 and M-6 Standard Missile [after 15] 

This chapter provides a brief overview of the 3-dimensional (3D) target model 
that operates in the Earth’s gravitational field, using available TPD-2 data from the open 
source. A 3-D interceptor model, based on the SM-6 specification, was developed that 
operates in the Earth’s gravitational field, was also used. Brief description of the various 
function files, assumptions made, data and values used in the simulation study is also 
presented to help the reader understand the program flow. A listing of the MATLAB 
program code as well as detailed description and theory of the various functions are 
reproduced in Appendices D and E, respectively, for completeness and easy cross- 
reference. 
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B. 


SIMULATION SOFTWARE ARCHITECTURE 


The general program flow of the 3DoF simulation conducted is shown in Figure 
7. The remaining sections of this chapter will briefly describe the function of the main 
blocks shown in the architecture, models used and the general flow of the simulation. 
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Figure 7. Comparison of the size of TPD-2 and M-6 Standard Missile [15] 

C. TARGET MODELING 

1. Ballistic Missile Physical Specifications 

Table 1 shows the physical characteristics and specification of the TPD-2 missile 


[13]. 



Overall 

Length 

32 m 

Payload 

750-1000 kg 

Range 

3500-4300 km 

Stages 

2 



Stage 1 

Stage 2 

Diameter 

2.2 m 

1.335 m 

Length 

16 m 

14 m 

Launch Weight 

-60,000 kg 

15,200 kg 

Thrust 

-103,000 kg f 

13,350 kgf 
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Fuel / Oxidizer 

TM-185/AK- 

271 

TM-185/AK- 

271 

Propellant Mass 

- 

12,912 kg 

Bum Time 

-125 s 

110s 


Thrust 

Chambers 

4,1 

Type 

LRICBM 


Table 1. Specifications of TPD-2 Ballistic Missile 

2. The Ballistic Missile Model Program 

The 3DoF model developed by Lukacs and Prof Yakimenko in 2006 comprises a 
series of MATLAB functions on an iterative integration loop, using 4 function files to 
accomplish the modeling. A brief description of each function is as follows: 

1. BRFlight3.m - integrates each time step to determine the current position, 
attitude, and aerodynamic forces acting on the rocket/missile; 

2. BRParams3.m - determines the mass of the ballistic missile and the 
surface reference area; 

3. BRDrag.m - determines the drag coefficient; 

4. STatmos.m - determines the properties of the local atmosphere; 


The program BRFlight3.m generates a ballistic flight path of the ballistic missile 
target based on the model developed by Zarchan [16]. The program generates a 3-D 
flight path that is contained within the 2-D x-z plane shown in Figure 7, where the 
asterisks represent the staging events. 



Y position (m) X position (m) 


Ballistic Missile Flight Path [Ref ] 
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Figure 8. 






















The launch position and angle are as follows: 

*0 =° 

To =0 
z 0 = Re 

<9 = 85° 

where Re is the radius of the Earth (6,378,137 m) based on WGS-84 system. 

The thrust is assumed fixed for each stage of the propulsion. At 130 and 240 
seconds after launch, the thrust drops instantaneous to represent the staging events. At 
completion of the boost phase, the thrust is zero. The thrust profile of the two stages is 
shown in Figure 8. 
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Figure 9. TPD-2 Thrust Profile (Boost Phase) [Ref] 

In order to account for the aerodynamic forces, drag has to be calculated. 
Atmospheric conditions, like density and temperature, are first determined using the 
function STatmos.m. The missile drag coefficient, CD is then computed by the function 
BRDragC.m. This is followed by calculating the drag force. 

The ballistic missile’s mass is a simple function of time (as shown in Figure 9) 
and this is computed using the function BRParams3.m. The mass drops sharply during 
the staging events at 130 and 240 seconds. After the completion of the boost phase, the 
mass remains constant for the remaining duration of the flight. 
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Figure 10. TPD-2 Ballistic missile Mass (Boost Phase Only) 

The final step of the ballistic missile model is to record all the data for the ballistic 
missile. This data will be called by the interceptor model to simulate detection by the 
missile’s onboard sensors. The interceptor missile will then register the location and 
velocity of the target missile at the appropriate intervals. In this program, a 60-second 
delay is assumed for the launch of the interceptor missile to account for the detection and 
decision loop. 

The final ballistic fly-out range and the altitude profile predicted by the model is 
presented in Figure 10 and 11. 



Figure 11. TPD-2 Fly-out Range 
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Figure 12. TPD-2 Altitude Profile for Entire Flight (left) and Boost Phase (right) 

D. INTERCEPTOR MODELING 

1. Basic Definitions and Assumptions 

The basic specifications of the SM-6 missile are presented in Table 2 [17]. 



Overall 

Length 

6.5 m 

Payload 

115 kg 

Range 

150 km 

Stages 

2 

Thrust 

Chambers 

1,1 

Type 

ERAAW 



Stage 1 

Stage 2 

Diameter 

0.53 m 

0.34 m 

Length 

1.72 m 

4.78 m 

Launch Weight 

712 kg 

686 kg 

Thrust 

- 

- 

Fuel / Oxidizer 

HTPB-AP 

TP-H1205/6 

Propellant Mass 

468 kg 

360 kg 

Bum Time 

6 s 

- 


Table 2. Specifications of SM-6 Interceptor Missile 

2. Interceptor Missile Model Program 

The 3DOF model used in the program comprises 4 MATLAB functions files ran 
on a repeating integration loop. A brief description of the function files are as follows: 

1. SMFlight3.m - integrates each time step to determine the current 
position, attitude, and aerodynamic forces acting on the missile; 
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2. SMParams3.m - determines the mass of the interceptor missile and 
the surface reference area; 

3. SMDrag.m - determines the drag coefficient; 

4. STatmos.m - determines the properties of the local atmosphere. 

Drag is calculated in a similar fashion as that implemented in the ballistic missile 
model. The SMParams.m function uses the same methodology as the BRParams.m 
function to calculate the reference surface area and mass. 

The thrust profile for the interceptor missile is also a two-step curve, just like the 
ballistic missile profile. The only difference is the timing for the staging event. This is 
shown in Figure 11. The staging event occurs at 6 and 26 seconds. 

2.5 


H 1 

0.5 


Time (sec) 

Figure 13. SM-6 Interceptor Thrust Profile 

The mass of the interceptor mission changes at 6 and 26 seconds to represent the 
staging events. Figure 13 shows the mass profile. It shows a distinct drop at 6s and a 
constant mass after 26s. 
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Figure 14. SM-6 Interceptor Mass (Boost Phase Only) 

E. INITIALIZATION 

For initialization, the earth geographical data required was based on the WGS84 
values presented in Table 3. [18]: 


Earth's Radius, R e 

6,378,137 m 

Earth's Semi-Minor Axis, b 

6,356,752 m 

Earth's Flattening (1/Ellipticity),/ 

1/298.257223563 

Earth's Rotation Rate (Qz E ) 

7.292116 e-5 rad/sec 

Earth's Gravitational Constant, GM 

3.986004418e 14 m 3 /s 2 


Table 3. WGS-84 Values 

F. COMMON FUNCTIONS 

Two functions were commonly by all the models. They are: 

1. SMDragC.m and BRDrag.m 

The drag on the missile is dependent on two conditions, the Mach number and the 
phase the missile—boost or glide. Mach number represents the speed of missile and the 
phase will determine the presence or absence of the base drag. The plot of drag 
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coefficient with Mach number of the missiles is shown in Figure 14 [19]. The same 
function is applied through two drag files in the interceptor and ballistic missile model 
respectively. 


Drag Coefficient by Mach Number and Boost Phase 



Figure 15. Drag Coefficient by Mach Number and Flight Phase 

Having obtained the drag coefficient and computed the dynamic pressure 
computed, drag can be calculated. 

2. STatmos.m 

The atmospheric parameters, which are dependent on altitude, are computed by 
STatmos.m function file. The calculation was based on the 1976 standard atmospheric 
survey, and includes values up to 86 km in a tabular format. Figures 15-17 show the 
atmospheric charts used by STatmos.m to derive all the parameters required. 


Atmospheric Temperature Variation by Altitude 
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Figure 16. 


Atmospheric Temperature Variation by Altitude 


Atmospheric Density Variation by Altitude 



Figure 17. Atmospheric Density by Altitude 


Atmospheric Pressure Variation by Altitude 



Figure 18. Atmospheric Pressure Variation by Altitude 

G. THE FLIGHT PROGRAM 

Once the previous iteration (or the initialization) has been integrated, it is returned 
to the program as the current values. The program then uses these values to calculate all 
the descriptive values of the system, apply the corrective time- and position- dependant 
factors, and calculate the derivatives for the next iteration. 

The BRParams.m and SMParams.m functions were called to determine the 
reference areas for the control surfaces and missile plan form areas. The program then 
calls the function ACoeff.m to determine several necessary aerodynamic coefficients. 
The function inputs are angle of attack, altitude, and pitch control surface deflection. 
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The respective missile model will call the function SMDrag.m and BRDrag.m to 
determine the value of the drag coefficient. Drag and thrust are then calculated and the 
execution is returned to the main program. 

The trajectory of the interceptor missile is generated by the function 
SMTrajectory.m by applying the TS guidance law through the SMGuidance.m function. 
This is an iterative process, starting with a ‘guess’ of the interceptor final states and 
subsequently performing trajectory optimization by calling sub-routine 
SMGuidanceCost.m to minimize the cost and penalty of each iterative flight path 
generated. The optimized flight path is then returned to SMTrajectory for execution 
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III. COMPARISON OF COMPENSATED PROPORTIONAL 
NAVIGATION WITH TRAJECTORY-SHAPING GUIDANCE 

In order to understand the guidance laws, it is useful to understand what 
“guidance” is. Guidance is the logic that issues steering commands to the vehicle to 
accomplish certain flight objectives’ [20], In the case of a missile and its target, it is the 
algorithm the missile uses to guide itself in an intercept path with the target, through 
autopilot commands given to its control surfaces, like fins or wings. In current interceptor 
missiles being fielded or deployed, a number of guidance laws were used. The type of 
guidance law implemented depends on the mission and the type of targets the missile is 
designed for. Examples of better-known guidance laws are the Beam Rider, Pure Pursuit, 
Proportional Navigation and its variants. 

A. PROPORTIONAL NAVIGATION GUIDANCE 

PN guidance or its variants, like augmented or compensated PN, are the most 
common guidance law implemented in modem missiles. PN guidance has proven itself in 
many operations to be effective against maneuvering target, especially in an air-to-air 
duel or surface-to-air interception. It can be found in many air-to-air, air-to-surface and 
surface-to-air applications. In proportional navigation, the fundamental principle lies in 
the rate of change of the missile heading being kept proportional to the rate of rotation of 
the line-of-sight (LOS) from the missile to the target. The guidance system points the 
differential velocity vector at the target [7]. Figure 18 illustrates the geometry of 
interception using the PN guidance. 
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Figure 19. Intercept Geometry for PN Guidance [after 21] 

In order for interception to occur, the missile heading angle, 'P and the LOS angle, 
X, must be kept constant. If the target increases speed, then X will increase and the 
interceptor must correspondingly increase ¥ to maintain its interception course. The rate 
of change of both angles must be proportional, given by the expression [16], 

't' = NA (III.A.l) 

where N is the proportionality constant, which depends on the target lead factor of the 
interceptor (usually between the values of 3-5 depending on the interceptor missile’s 
maneuverability). The interceptor will maneuver until A = 0 (no more changes in LOS 
rate with the target), at which point 'P = 0 . The interceptor continuously adjusts itself to 
maintain the above condition till interception (or a miss) occurs. 

During the intercept process, the missile seeker continuously points itself at the 
target to derive the LOS rate and provides a signal to the guidance computer. The 
guidance computer in turn will convert the feedback into an acceleration command to the 
control surfaces to physically steer the missile. The normal acceleration to the missile’s 
velocity vector is given by 

a = V Dx ¥ (III.A.2) 

which can be correlated to the LOS rate 

a = NV°A (III.A.3) 
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This is applicable to a 2-D case. Zipfel [20] develops a 3-D dimensional PN guidance law 
with the expression 

a = NV d Q oe u v - g (III.A.4) 

where A is the cross product of the LOS frame with respect to the inertia Earth frame, 
u v is the unit vector of V D , and g is the added gravity bias, which counteracts the sagging 

tendency of the trajectory under seeker control. This form of guidance is termed by Zipfel 
as Pure PN. 

B. COMPENSATED PROPORTIONAL NAVIGATION 

Pure PN, however, is not commonly implemented in its original form. The thrust 
generated by the interceptor’s propulsion creates a parasitic acceleration in the LOS angle 
that has to be compensated in the autopilot command. Guidance in this form is termed 
compensated PN. The parasitic acceleration projected onto the LOS plane can be 
expressed as follows: 

[a m r = WS M (III.B.1) 

It is then subtracted from the PN command in equation (3.B.4) to obtain the compensated 
command [20] 

a = NV D [£2«]'[<,.]" - ws w R[a m ]" - ~,R [g]' (IH B.2) 

where the rotation matrix L ™R of the LOS coordinates with respect to the wind frame is 
defined by the azimuth and elevation angles from the LOS vector. 

C. DISADVANTAGES OF COMPENSATED PN GUIDANCE FOR BOOST 

PHASE INTERCEPTION 

Three disadvantages of the compensated PN for the boost phase interception of 
ballistic missiles were identified by Lukacs and Prof Yakimenko. The first disadvantage 
is the inherent control system time constant of the PN guidance, which utilizes current 
target information in a homing guidance loop with feedback control. Such control loops 
will inevitably result in a finite time constant which can result in considerable miss 
distance. The second is the disregard of the end-game environment. The PN guidance law 
is focused on maintaining the LOS rate constant continuously based on current missile 
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and target parameters. However, the guidance law does not attempt to optimize missile 
speed and altitude, which is critical for the endgame condition, to ensure that the missile 
interceptor has sufficient acceleration left to intercept the target. Rather, PN guidance 
adopts the most direct and minimum acceleration collision path with the hope that there is 
little target maneuver and hence sufficient acceleration left in the interceptor for the 
endgame. Lastly, the intercept geometry is not controlled for the case of PN guidance, 
which is left to the relative speeds of the interceptor and target as well as the endgame 
maneuvers involved Controlling the intercept geometry allows for the maximization of 
kinetic energy transfer to the target, which is critical for boost phase intercept of ballistic 
missile since a lightweight (very small or no warhead), high speed missile with sufficient 
range (most of the weight goes to the fuel) is required. 

D. ADVANTAGES OF TRAJECTORY SHAPING GUIDANCE 

In TS guidance, the derivation of the intercept solution - the optimized interceptor 
flight path, results in a set of functions, which occurs through a Cost Function and a 
Penalty function, which must be minimized through multiple iterations. The three 
parameters represented in the Cost Function are the length of the virtual arc (proportional 
to the flight path distance), the time to intercept and the impact angle of the final intercept 
(angular deviation from desired 90° impact). Each of the parameters needs to be carefully 
weighted to correctly reflect the desired condition of the intercept, and affects the Cost 
Function value as a whole. 

The penalty function consists of the maximum acceleration in the y- and z- 
directions. They represent the ‘penalty’ to pay when certain physical limitations are 
reached or violated. The penalty function has been set to the certain values, which are 
dependent on speed and altitude of the interceptor, to represent the physical limits beyond 
which the intercept solution will not be feasible. 

The flight path optimization is done by solving a non-linear programming 
problem numerically real-time and once the minimum function is obtained, the algorithm 
will return the required control time history to the missile guidance system, which can 
then execute the commands and fly the derived flight path. Since the missile system can 
be programmed with sufficient data to compensate for its control system time constant, 
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the system lag can be effectively negated, thus eliminating a source of error. The 
guidance system can be updated every few seconds to increase the accuracy of the 
intercept. 

The main advantages of this guidance are three-fold—(1) the relatively short 
computation time to iterate and converge to an optimized solution, (2) the cost and 
penalty functions are scalable as desired to fit the mission profile and (3) elimination of 
the control system time constant. For a boost phase intercept mission, the TS guidance is 
able to address the disadvantages of the computed PN guidance and offers greater 
flexibility in being able to ‘customize’ the guidance to improve its performance to meet 
the different operational demands. 
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IV. INDUCED DRAG MODEL 


A. AERODYNAMIC DRAG 

The modeling of aerodynamic forces acting on a missile is a necessary step in the 
simulation of the missile interception of an aerial target. Drag, which represents the total 
sum of forces that opposes the thrust provided by the missile propulsion system, is a 
‘major design parameter in satisfying the flight range requirement of tactical missiles’ 
[22], This is especially so for supersonic missile. It is part of the guidance algorithm to 
incorporate a drag model to account for the effect of drag on the missile, which affect 
flight performance in all phases. A suitable drag model will enhance the accuracy of the 
simulation result and help to provide a more realistic prediction of the flight characteristic 
and performance. 

B. COMPONENTS OF DRAG 

In the earth atmosphere, any flying object will experience drag. It is the force that 
opposes thrust (or the forward motion) of aircraft, missiles or other flying platforms. The 
overall drag on a body is made up of several components—friction, base, wave (or 
pressure) and induced drag. Friction drag refers to the drag generated by the skin friction 
between tangential air flow and the moving body. It is primarily a function of the body 
fineness ratio or length-to-diameter ratio of the missile body. Base drag accounts for the 
pressure difference of air flow between the nose and the base of the body. It is 
particularly significant during coast or glide phase of a missile since the propulsion 
system is no longer generating thrust. Wave drag is caused by the effect of normal 
pressure of air and largely dependent on the shape of the nose. The nose fineness ratio 
gives a measure of the nose shape, which affects the wave drag. A sharp nose will create 
less drag than a blunt one. Last but not least is induced drag. This is the component of 
drag that is ‘induced’ by the lift generated for level flight or climb. It is directly 
proportional to the lift. Hence induced drag is present in all flying platform as long as lift 
is generated by the control surfaces and in some condition, is a key contributor to the 
overall drag. 
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In mathematical terms, drag is a function of the drag coefficient, dynamic 
pressure and reference area, given by this equation 


(IV.B.l) 


D m 


The drag coefficient, Cd is usually derived from empirical equation as a function 
of the speed of sound or Mach Number (M) in the operating regime, the angle of attack, 
a, which is the angle of inclination between the missile body axis and its velocity vector. 
In the case of a missile, the reference area is commonly taken as the cross-sectional area. 

The atmospheric dynamic pressure, q is given by 



(IV.B.l) 


where p is the density of air and V is the speed of the missile. 

C. INDUCED DRAG POLAR 

As briefly described in Chapter 2, for the new TS guidance, a standard drag model 
was used whereby the equation (IV.B.l) was used to calculate the drag. Noticed it was 
mentioned that Cd is usually obtained empirically and different missile design is likely to 
have different drag characteristics. The accuracy of the drag prediction can be enhanced 
with the use of an induced drag model. This model predicts drag by involving the lift 
coefficient, Cl which describes the lift generated by the missile. The lift force, L is 
normal to the missile velocity vector and hence also normal to drag, D (which is opposite 
to the velocity vector). The lift and drag coefficient can be expressed as 



& 


(IV.C.l) 


And 



P 

G-ref 


(IV.C.2) 


Both coefficients can be assumed to be functions of the following parameters: 
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Cl, Cd =/ (Mach number, angel of attack, power on/off, shape) [20] 

The Mach number can have a great effect on the coefficients, especially in the 
transonic and supersonic regime. Of note, at the transonic regime, there is a sudden steep 
rise in drag before the speed crosses the sonic line, a phenomenon commonly known as 
transonic drag rise (shown in Figure 19). The main effect, however, is caused by the 
angle of attack. A slight change can result in significant change in lift coefficient. 



Figure 20. Transonic Drag Rise 


When the lift and drag coefficient are plotted against each other, a near parabolic 
curve emerges, for any given Mach number. The relationship can be expressed as: 

C D = C D0 + k(C L -C L0 ) 2 (IV.C.3) 

Graphically, the plot is shown in Figure 19 and is known as the parabolic drag 

polar: 


C L 



Figure 21. Parabolic Drag Polar for Asymmetric Body 
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In the drag polar expression, Cdo is the zero-lift drag coefficient, a term that 
accounts for the other drag components when lift is zero while Clo represents the lift 
coefficient at the lowest drag. The coefficient k is referred to as the induced drag 
coefficient. 

If minimum drag occurs at zero lift, then Clo is zero and the parabola is 
symmetrical about the Cd axis and centered at Clo = 0. This case is representative of a 
missile, which is axis-symmetric in nature. The different points on the drag polar 
represent the different angles of attack, which give rise to different lift and drag 
coefficient. The parabolic drag polar of a missile is as shown in Figure 20: 



Figure 22. Parabolic Drag Polar for Symmetric Body 

In order to implement the induced drag model in the algorithm to calculate drag 
(D), the set of Cl-Cd data for the interceptor and target at various Mach number was 
obtained from open sources and set up in the form of a look-up table. The data is then 
curve fitted using a polynomial (parabolic) function. A plot of the family of curves for 
different Mach number using Excel plot is shown in Figure 21: 
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Figure 23. Parabolic Drag Polar for Different Mach Number (Symmetric Body) 


The values of Cdo and k can be derived from the coefficient of the polynomial 
function. These values can be substituted into equation (IV.C.3). Since Clo is zero for an 
axis-symmetric missile, the only outstanding variable that has to be computed is the lift 
coefficient, Cl. For a particular iteration, Cl can be expressed as a function of load factor, 
n z . The load factor of a missile is defined as: 


n z 


L_ 

w 


(IV.C.4) 


Where L is the lift and W is the weight of the missile. From equation (1V.C.1), we 
can express L in terms of Cl. 

L = C L qS re f (1V.C.5) 

Therefore, substituting equations (1V.C.5) into (1V.C.4), Cl can be computed: 




(1V.C.6) 


The overall drag coefficient can then be calculated using equation (1V.C.3), with 
values of Cdo, k and Cl. 

A MATLAB program was developed to implement the induced drag model in the 
TS guidance algorithm. The code is appended in Appendix C. Minor modifications have 
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to be made in the other sub-routines in order to integrate the new drag model into the 
algorithm. It is intended to replace the generic drag model currently. Though the new 
model has been integrated and initially tested to be working, results has not been 
consistent. Further testing using the new code and more simulation scenario to fine-tune 
the program can be conducted in future studies. 
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V. RESULTS AND DISCUSSIONS 


A total of 270 simulations were conducted in this trade-off study. The main 
objectives of the simulation study are: (1) to verify that the TS guidance can be used in a 
hit-to-kill interceptor missile to effectively engage a ballistic missile in the boost phase, 
(2) to predict the region of the interceptor launch position in relations to the ballistic 
missile launch point which can result in a feasible intercept solution and (3) the impact on 
the time-to-intercept, altitude and impact angle deviation by varying the intercept 
geometry coefficient, Wia of the cost function. In this study, simulations were run for 
various interceptor missile (simulated by the surface-launch SM-6 missile) launch 
positions, for different northing and easting against a ballistic missile (simulated by the 
TPD-2 ICBM) launched from a fixed site (arbitrarily fixed at point [0,0,Re] in the earth 
coordinate system, whereby Re is the radius of the earth). The ballistic missile is 
launched towards a target in the north direction. 

A. FEASIBILITY OF THE TRAJECTORY SHAPING GUIDANCE 

The simulation results showed that the trajectory guidance algorithm worked well 
for the intercept scenario ran and are within expectation. The guidance provides feasible 
solutions for the interception of the TPD-2 ballistic missile by the surface launched SM-6 
missile. Out of the 270 simulation runs, successful intercept was achieved for 180 runs. 
The remaining ‘non-feasible’ intercepts were due to the inability of the guidance 
algorithm to converge to an optimized problem within the maximum number of iterations 
allowed. The maximum iterations condition was set to act as the upper bound and ‘time¬ 
out’ condition, reflection the real world situation whereby the intercept scenario does not 
offer a solution fast enough to ensure a successful intercept. This is realistic as the 
interceptor missile has to compute the optimized flight path after launch and be guided to 
the predicted intercept point within a very small time window. Such scenario occurs 
when the interceptor system is ‘out of range’ or within the ‘blind range’ to effectively 
intercept the ballistic missile within the timeframe of 60 to 180 seconds after launch. 
Table 22 to 24 presents a summary of simulation results, showing the number of 
iterations the algorithm takes to find the optimized intercept flight path for values of Wi A 


33 



set to 0, 100 and 1000 respectively. The region in green represents feasible intercepts 
predicted by the TS algorithm and the numbers refers to the number of iterations required 
to find the optimized intercept flight path. The letter ‘EX’ in the red region represented 
the non-feasible region whereby there is no feasible intercept solution if the interceptor 
missile is launched from those positions. 



Table 4. No of Iterations Required for Feasible Interceptor Solution (WIA = 0) 



Table 5. No of Iterations Required for Feasible Interceptor Solution (WIA = 100) 
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Table 6. No of Iterations Required for Feasible Interceptor Solution (Wi A = 1000) 


The results obtained successfully verified that the TS guidance algorithm does 
provide feasible solution to the boost phase intercept problem. The new guidance is 
suitable for scenario involving a high speed hit-to-kill missile interceptor and a fixed 
trajectory ballistic missile target. By comparing the results for the different intercept 
geometry coefficient values, it is observed that the number of iterations is increased 
across the board when Wia is increased from 0 to 1000. This means that when a higher 
weightage is placed on intercept geometry in the intercept solution, it increases the 
computing resources required and the optimization process takes a longer time to 
converge. 

B. REGION OF POSSIBLE DEPLOYMENT POSITION OF THE 
INTERCEPTOR MISSILE 

The simulation study was able to provide a mapping of the feasible intercept 
region, in terms of the interceptor missile launch position with respect to that of the 
ballistic missile. As the scenario is symmetrical about the north-south direction, only half 
of the 360-degree coverage needs to be tested. The other half is identical and a mirror 
image of the results obtained. The simulation results obtained and plot is appended in 
Appendix 1. Figure 22 shows the feasible intercept region in green and non-feasible 
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region in red, for the scenario tested. The mapping was done by running simulations and 
collecting intercept data at intervals of 20 km in the northing and easting direction. 



- 140 - 120-100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 


Easting (km) 


Figure 24. Region of Interceptor Position for Interception of Ballistic Missile 

It should be noted that the shape of the boundary has no special significance 
except for the fact it is purely a 2-D plotting limitation of the software used. A smooth 
curve should be the correct representation for the boundary of the respective regions. 
From the mapping, several interesting observations were noted. As expected, the plot 
shows an off-set in the feasible intercept region with respect to the ballistic missile launch 
point. This represents a distinct difference in the intercept range between front- or rear- 
quarter engagement scenarios. One would expect the rear-sector engagement, represented 
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by the intercept launch points south of ballistic missile launch point at [0,0], to have a 
feasible region with shorter ranges. This is consistent with a typical ‘tail-chase’ intercept 
scenario, which requires the interceptor to be closer to the target. The front sector 
engagement is characterized by longer ranges as depicted by the larger area in the plot. 

There is a central non-feasible region (in red) and for near ranges close to the 
ballistic missile launch point. This indicated that certain aerodynamic limits were 
exceeded, like the g- and lateral acceleration limits, which resulted in a high penalty 
function value and divergence during the optimization process. The latter represents some 
kind of ‘blind-range’ whereby intercept is not feasible due to the target being too close 
for the interceptor to effectively maneuver and to achieve an intercept. 

The usefulness of mapping the region of feasible intercept is that it provides 
critical information for the operational planners on the possible region to deploy an 
interceptor system (using the TS guidance) given a ballistic missile threat. The creation of 
such plots can be automated and done in a short time using computer program. The input 
requirements are the specifications of the interceptor missile and intelligence on the 
specifications/characteristics of the target missile. The accuracy of such prediction 
software would depend on the number of simulation run (proportional to the ‘resolution’ 
or distance interval required) and accuracy of the input data and model used. There is 
scope to develop a computer program to generate intercept region plots for operational 
planning purposes in future studies. 

C. EFFECT OF VARYING THE INTERCEPT GEOMETRY COEFFICIENT 

One of the main objectives of the trajectory-shaping guidance is to effect 
maximum transfer of kinetic energy from the interceptor missile to the target missile in a 
hit-to-kill interception scenario. This can be achieved when the impact angle onto the 
target missile is at a right angle to the ballistic missile body. The TS guidance algorithm 
allows the missile designer to adjust the weightage placed on the intercept geometry, 
amongst other critical parameters, by varying the value of Wia, which is to be optimized 
(minimized). The trade-off study examines the effect on time-to-go, intercept altitude and 
impact angle deviation due to changes in Wia. 
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In this study, simulations were run for three values of Wia, namely 0, 100 and 
1000. By having Wia= 0 implies that the optimized flight path does not take into account 
the end-game intercept geometry. In the case of Wia = 1000, a higher weightage is placed 
on impact angle at intercept and the flight path will be optimized to achieve an impact 
angle as close to 90 degrees as possible. When a higher ‘premium’ is placed on the 
intercept geometry, the algorithm will attempt it at the expense of other performance 
parameter. Here, there is always some form of trade-off that a missile designer has to be 
conscious of. 

A comparison of the time-to-go or time-to-intercept, intercept altitude and impact 
angle deviation for different values of Wia provides valuable insights into the effects of 
Wia on the missile performance under the same intercept scenario. Table 7 presents a 
summary of the comparison results. 







Time-to-intercept (s) 

Intercept Altitude (km) 

Impact Angle Dev (°) 

WIA 

Min 

Max 

Ave 

Min 

Max 

Ave 

Min 

Max 

Ave 

0 

17.6 

61.0 

38.67 

22.3 

60.0 

38.77 

2.4 

30.3 

18.45 

100 

17.6 

61.0 

38.48 

22.3 

60.0 

38.77 

2.4 

30.3 

18.45 

1000 

17.8 

60.2 

39.00 

22.4 

59.1 

38.94 

0.0 

30.0 

10.02 


Table 7. Time-To-Intercept, Altitude And Impact Angle Deviation 
For Different Intercept Geometry Coefficient, Wia 

Generally, it can be seen that there is no significant different between the results 
obtained for Wia being set at 0 and 100, for the time-to-intercept, intercept altitude and 
impact angle deviation. However, some different is observed when the weightage 
coefficient is increased by a factor of 1000. The most significant difference is observed in 
the impact angle deviation. With a high weightage placed on intercept geometry (Wia = 
1000), the minimum impact angle deviation is 0, and the average over all the feasible 
intercepts reduces from about 18° to 10°' This is within expectation and verified that the 
algorithm indeed attempted to minimize the impact angle deviation when the Wia is 
increased. Another important observation is that by enhancing the intercept geometry 
weightage, it does not affect time-to-intercept and intercept altitude. By correlating the 
observations made in the previous section, it can be seen that when Wia is increased, the 


38 




number of iterations needed for convergence also goes up. However, the overall time 
taken for the intercept generally does not increase significantly. This seems possible if 
there is some trade-off or interaction with the other parameter(s) which is not investigated 
in this study. It can also be seen that the area of intercept decreases when Wia is changed 
from 0 to 1000. A possible reason is that the larger number of iterations required for 
higher intercept geometry coefficient may have resulted in the further range intercepts 
exceeding the maximum iteration limit imposed by the guidance. This will cause the 
longer range region to be marked as non-feasible. However, this clearly shows the 
‘penalty’ for imposing a placing a higher weightage on intercept geometry. Figure 23 
provides an illustration of the typical intercept geometry at Wia = 1000. It can be seen 
that the flight path is optimized for a near 90° impact to the ballistic missile trajectory. 



Figure 25. Typical Intercept Geometry 

It can be expected that when the weightage for intercept geometry is increased, 
there is greater constraints imposed on the interceptor system to optimized impact angle 
amongst other factors. This will affect missile performance in terms of time-to-go and 
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higher intercept altitude. A quick sensitivity study by using WIA values of 0, 100 and 
1000 and computing the statistical average of all feasible intercept for time-to-go, 
intercept altitude and impact angle deviation in the study revealed that increasing the 
weightage by a factor of 100 does not result in significant changes in the missile 
performance. At a factor of 1000, the effect on the missile flight parameters begins to be 
more observable. More simulations can be conducted in future studies for higher values 
of the Wia that will cause significant effect on missile performance. Some form of scaling 
factor can be obtained, from which a missile designer can select suitable values Wia to 
use for different missions. 

There is hence a trade-off between the weightage coefficient and other 
performance parameters of the interceptor missile. From the trade-off study, it was 
demonstrated that the TS guidance allows the missile designer the flexibility to adjust the 
weightage of the critical parameters in order to optimize the intercept path. With the 
insight obtained from the study, the effect of the weightage coefficient on missile 
performance is known. There is a need to balance the trade-off between the weightage 
placed on the different critical parameters. In the particular case of intercept geometry, 
the weightage coefficient (Wia) must be carefully selected to optimize missile 
performance for a particular mission or scenario. The same consideration can be 
translated to the other critical parameters. Again, a computer program can be developed 
to aid missile designer in selecting a suitable weightage or cost function coefficients for 
different mission requirement and operational scenario. This will be very useful if TS 
guidance in implemented in intercept missiles and offers greater advantage over other 
guidance laws as one that allows ‘customization’ for different mission requirements 
simply using software changes. 


D. INTERCEPT GEOMETRY TRENDS 

The simulation study also provides good data to perform simple trending analysis. Useful 
charts can be generated to provide trending information. An example is the impact angle 
deviation data derived from the simulation. It provides greater insights into how the 
intercept geometry will change with respect the ballistic missile trajectory azimuth angle, 
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0bm, with respect to the interceptor launch site. The angle 0bm, is defined as the angle 
between the trajectory plane of the ballistic missile and the direction from the interceptor 
launch point to the detected ballistic missile position, as shown in Figure 26. 



Traiectorv 


A 



Ballistic 
Missile Launch 
Direction 


Ballistic 
Missile 
Launch Point 


Figure 26. Illustration of Angle between interceptor Launch Position and 

Ballistic Missile Trajectory Plane 

Figure 27 shows the plot of impact angle deviation for interceptors launched from 
the various positions. 





Easting, km 


Figure 27. Effect of Interceptor Launch Direction on Intercept Geometry 
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It can be observed from the plot that interceptors launched normal to the ballistic 
missile trajectory plane (0bm = 90°) usually have small impact angle deviation (very close 
to 0°) or impact the ballistic missile at close to 90°' The impact angle deviation increases 
when the interceptor flight path becomes more oblique to the ballistic missile trajectory 
plane (that is, 0 B m becomes smaller). For example, the impact angle deviation is larger, 
when the interceptor is fired from [80 km N, 0 km E], as compared to an interceptor 
missile being launched from [0 km N, -50km E] 

The largest impact angle deviation occurs when the interceptor missile is in the 
front-sector of the ballistic missile. This can be explained from the fact that the missile 
does not have sufficient time to effectively ‘shape’ its flight path for a 90o impact on the 
ballistic missile, compared to the cases where the interceptor is closing-in from the side 
(normal to the ballistic missile trajectory plane). By plotting the results for the other 
parameters and analyzing them in future studies, other useful insight into the interceptor 
missile performance can be deduced. This will be helpful to missile designers as well as 
operational planners to optimize the potential of the TS guidance for possible anti- 
ballistic missile defense missions. 
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VI. CONCLUSIONS 


The simulation study verified that the TS guidance provides feasible solutions to 
the boost phase intercept problem involving the interception of a fixed-trajectory ballistic 
missile by a surface-launched interceptor missile. For the particular scenario used in the 
study, it was shown that if an interceptor system using the TS guidance can be deployed 
within 30 to 60 km (for front sector intercept) and 20 - 50 km for (rear-sector intercept), 
it is possible to intercept the ballistic missile effectively during the boost phase. The 
results also shows that the new guidance allows missile designer to ‘customize’ the 
guidance algorithm for specific mission by changing the ‘weights’ of some critical 
parameters However, there is a trade-off between the each of the weightage coefficient 
and missile performance. By enhancing the weightage on intercept geometry, there is no 
change in the time-to-intercept or intercept altitude. However, the region for the feasible 
interceptor launch position around the ballistic missile launch point has reduced. Hence, 
the weightage coefficient will have to be carefully chosen to ensure that there is a balance 
of benefit and penalty with respect to specific missions. This aspect of ‘customization’ 
represents a clear advantage of the TS guidance over other guidance law. The new 
guidance algorithm can be further enhanced and fine-tuned. Further research can include 
conducting simulation study using a 6 DoF model of the TS guidance, combining 
guidance solution with on-board navigation solution, testing the complete guidance, 
navigation and control solutions and developing a computer program as a design tool to 
perform trade-off study and select weightage coefficients for different mission 
requirements. The new TS guidance holds promise to be a feasible and implementable 
solution to the boost phase ballistic missile intercept problem that will become more 
prominent in the coming years. 
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APPENDIX A: TABLES AND PLOTS OF SIMULATION RESULTS 
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APPENDIX B: PARTICIPATION IN AIAA YOUNG 
PROFESSIONAL ROCKET LAUNCH COMPETITION IN MAY 2009 


INTRODUCTION 

The author participated in the AIAA Young Professional Rocket Launch 
Competition in May 2009 as a member of the NPS Team (see Figure 28 for the team 
photo). Team Peacock, one of the two teams from NPS, comprised eight student 
members and a faculty advisor (Prof Yakimenko). The competition required its 
participant to apply their knowledge to build a rocket, analyze and predict it performance 
and launch it at a test site, using common rocket kits, a choice of propulsion units and a 
payload provided by the organizers. The 8.5-feet long scaled model rocket was expected 
to reach height in excess of 7,000 at a speed of greater than Mach 1. The rocket has to be 
recovered after reaching its maximum height. 



Figure 28. Team Members and Advisor with Their Rocket at the Launch Site 

RELEVANCE 

The team’s participation in the competition is relevance to their course of study in 
NPS in general in that it provided an excellent opportunity for the students to apply their 
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knowledge in aerodynamics, propulsion, missile design and manufacturing in analyzing a 
problems and physically implementing the solutions to achieve an actual task or a 
mission—in this case building the rocket and launching it. In the process, greater 
experience and knowledge was gained, not only in the academic areas, but also in terms 
of project management, teamwork and communication. This is a good example of the 
‘hands-on’ approach towards learning and an interesting and impactful means of 
education. With respect to this project, this rocket is like a scaled model of an interceptor 
missile designed to intercept a ballistic missiles in the boost phase. The characteristics of 
the rocket resembles that of a hit-to-kill interceptor, vis high speed, Tight weight’ 
(includes only the body and fins, actuators, guidance and control), small payload (perhaps 
a small KE warhead or no warhead at all) and Tong range’ (high ratio of fuel mass to 
total mass). A 2-stage rocket would bear greater resemblance to a typical missile 
interceptor. The valuable experience gained provided the author greater insight and better 
understanding of the fundamental issues facing the design of interceptor missiles, 
particularly the aerodynamics, propulsion, stability and design trade-offs. 


COMPETITION 

The judging criteria for the competition are the maximum height reached by the 
rocket and the accuracy of the predicted performance versus actual performance. There is 
a winner each for the highest height reached and the closest predicted-to-actual 
performance category. Besides the actual launch itself, a report submitted by each team 
on their analysis and work done is also assessed to determine the winner for the latter. 


GOALS 

Team Peacock set itself two goals for this competition : (1) it aimed to be the 
winner for the highest altitude category and (2) to build the most stable rocket in flight 
and be able to recover the rocket successfully. This goal guided the team in its design and 
construction of the rocket. 
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TASKS 


The team took about 4 months to build the rocket and prepare for the launch, 
which was held in end May 2009. After the team was formed in Jan 09, it went straight to 
work by brainstorming what are the tasks to be completed, the timeline, resources 
required (including budget) and roles of individual members. Taking into account the 
goals and tasks identified, allocation of tasks and a schedule broad was drafted. 


The broad tasks facing the team from the start to launch were as follows: 

a. Budgeting . The team worked within a budget of $1000 to purchase all the 
required materials that is not provided by the organizer (only the rocket kit, propulsion 
unit and the common-use payload is provided). Table 8 provides a breakdown of the 
funds needed (for supplies). 


Item 

Location 

Cost($) 

Qty 

Total($) 

60” TAC-1 Parachute 

Giant Leap Rocketry 

86 

1 

86 

TAC-1 Kevlar Bag 

Giant Leap Rocketry 

25 

1 

25 

Tubular Kevlar Shock Cord (1/2" x 15') 

Giant Leap Rocketry 

25 

2 

50 

1/4" Eyelet Swivel 

Giant Leap Rocketry 

3 

4 

12 

AeroTech Electronic Fwd Enclosure 

Aerotech-Rocketry 

170 

1 

170 

Construction Materials 

Perfectflite.com 

100 

1 

100 

K270W(engine) 

Aerotech-Rocketry 

150 

2 

300 

Misc 




217 


TOTAL 960 


Table 8. Budget for Rocket Launch Competition 


b. Scheduling. An initial schedule was drafted (see Figure 30) so that members 
have a good sense of the key milestones and timeline for the various preparatory 
activities. The schedule was updated as project progresses to reflect the status of the 
various tasks. 
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Figure 29. Timeline for the Rocket Launch Project 


c. Resource Planning and Purchasing (Supply) . Some of the team members 

were assigned the role to plan and decide the items, the specifications and 
quantities needed as well as to source and purchase the items (the actual 
purchasing is done by the administrative staff) 

d. Structure Design . The rocket kit is provided by the organizers. It uses the 
Quasar 1/16 rocket, which comes with the nose cone, body sections, fins and also 
a parachute for recovery. The team used commercial software RockSim to predict 
the CG and CP calculation as well as to determine the structural dimensions for 
the rocket design. Figure 31 presents the CG and CP position and the static 
margin. Based on the calculation and prediction, the team decided to reduce the 
length of the rocket by 14 inches as well as halved the fin area to reduce the static 
margin to 2 and optimize the performance of the rocket so as to achieve the 
required stability and highest altitude. 
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Quasar Scale: 1/1 6 

Rocket length: 101.750 In., diameter: 4.QQ0 In., span diameter: 20.500 In. 
Rocket mass 165.SB? oz. , Selected stage mass 165.982 oz. 

Shown w/o Engines. 



Method CG In. CP In. CNa Static margin Analysis 

Barrowman 66.791 89.451 34.881 5.67 The rocket is over stable. 

RockSirn 66.791 90.307 42.143 5.88 The rocket is over stable. 

Figure 30. The Specifications and Dimensions of the Quasar Rocket Kit 

e. Rocket Motor Selection . Two engines, both with about the same impulse, 
were provided for selection. The team can choose either the engine that has higher 
thrust and a shorter bum time (K800) or the other with lower thrust but longer 
bum time (K270). After conducting simulation study of the predicted flight 
profile, based on the rocket modified dimensions, using RockSim and a 
MATLAB program written by one of the team member, it was found that the 
K270 best meet our needs. The engine motor specifications are as follows: 


Motor 

AeroTech K270 

Diameter (mm) 

54.0 

Length (cm) 

57.9 

Prop. Weight (g) 

1,188.0 

Total Weight (g) 

2,100. 

Avg. Thrust (N) 

247.6 

Max. Thrust (N) 

425.7 

Tot. Impulse (Ns) 

2,154.9 

Burn Time (s) 

8.7 


The motor performance chart for the K270 are shown in Figure 32. 
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Figure 31. The Specifications and Dimensions of the Quasar Rocket Kit 

f. Analysis (Simulation and Modeling). This is one of the key components of 

this rocket launch project. The team has to use various methods to analyze and 
predict the flight performance/profile of the rocket to be built. The team relied on 
several methods to analyze the flight performance of the rocket, including simple 
hand calculations for rough prediction, using RockSim, a commercial rocket 
software for amateur rocketeers and a student-developed MATLAB program 
(with Simulink model) written by one of our members. The analysis and 
simulation done helped the team to determine the structural modification needed, 
weight and balance, CG and CP position, stability, selection of engine and the 
maximum altitude prediction. Figure 33 shows the Simulink model. Comparison 
of results, derived from RockSim and own MATLAB model, can be made and it 
was found that they are within about 10% of each other. RockSim predicted a 
higher maximum height reached of 7,000 feet, while the Simulink model 
predicted 6,400 feet. 


62 


























Figure 32. The Specifications and Dimensions of the Quasar Rocket Kit 


g. Testing . Besides analysis, some of the rocket components have to be 
physically tested to have a higher assurance that it worked as designed and modifications 
made did not introduce problems to the design. The team tested the rocket engine at the 
propulsion laboratory to match the actual performance with manufacturer specifications, 
the Electronic Forward Closure (EFC) to test the trigger mechanism for the explosive 
charge to deploy the parachute and the actual parachute to select the one to be used for 
the rocket recovery system. During the first propulsion unit test, there was a bum-through 
of the engine casing due to a wrong part provided by the manufacturer. The problem was 
raised to the manufacturer and a replacement part was delivered. The discovery was an 
important one as it prevented the same occurrence from happening during the actual 
launch for our team as well as for other teams. The second firing was very successful and 
it verified the same propulsion characteristic given by the manufacturer. The EFC test 
verified that the trigger mechanism worked and would cause the separation of the rocket 
body after apogee to deploy the recovery chute. The parachute testing helped the team to 
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decide using a single chute system for recovery instead of a dual-chute system. Figure 34 
shows pictures of testing being conducted. 



Figure 33. Team Members and Advisor with Their Rocket at the Launch Site 

h. Construction . The actual construction of the rocket took about three 

weeks, after all the parts and materials were delivered. For a start, a rocket stand or jig 
was build to facilitate the assembling work. For this rocket, the body was made up of 
several cylindrical sections of thick cardboard material and the nose cone is made of 
plastic. The construction process begins with strengthening the body with layers of 
epoxy. The internal compartments and fins were joined to the body using fibre-glass cloth 
with epoxy and left to cure. The various structures were attached to the rocket body or 
inserted inside the rocket in sections. Fixing the fins to the tail section required more 
attention to ensure they are symmetrically attached. Much effort was put in to ensure the 
rigidity as the fins as they would be subjected to relatively large forces in flight. The 
rocket engine casing, EFC and parachute were then put inside the rocket body. Once all 
the components were assembled in their sections, the three sections were then joined 
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together using internal connectors and nose cone was attached to the body. The last step 
is painting the rocket to give it a distinctive color for aesthetic reasons as well as to 
provide a smooth finishing to reduce skin drag. The payload would be placed inside the 
rocket only on launch day. It consisted of an altimeter that would send altitude data 
continuously to the ground receiver for recording purposes. Figure 37 shows some 
pictures of the rocket being built at the laboratory. 



Figure 34. Team Members and Advisor with Their Rocket at the Launch Site 


i. Launch . The launch was the climax of the five-month effort and the team 
was as prepared as it could be. It was held on 26 May 09, Saturday at Koehn Lake 
Launch Site near to Mojave. On that day, there were 9 teams present with their rockets. 
Team Peacock was the first team to launch. After the rocket propellant was inserted into 
the casing, it was mounted onto the launch rail for the ignition. The igniter was fired by 
the firer, who was positioned at the firing ‘bunker’. All other observers were herded to 
the observation deck to watch the firing. At the countdown of 10, the rocket was 
launched. It lifted-off vertically with a sharp ‘bang’ and within seconds, the rocket motor 
accelerated the rocket to more than 3,000 ft and out of visual sight, leaving a trail of 
smoke behind it. Moments later, the parachute was deployed and the rocket descended to 
the ground in two parts, as expected. It landed ‘safely’ at the desert ground close to the 
launch site and was retrieved by the recovery team. Figure 36 shows the launch of the 
rocket. 
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Figure 35. Team Members and Advisor with Their Rocket at the Launch Site 

RESULTS 

The results made all the hard work pay off. Team Peacock emerged the winner of 
the closest altitude prediction category. The recorded maximum altitude reached during 
the actual launch was 5,700 feet. This was about 9% deviation from the predicted height 
based on the team’s analysis. Figure 37 presents the results in terms of the altitude versus 
time plot that was recorded by the organizer. 
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Figure 36. Official Altitude Plot for Team Peacock 


LESSON LEARNED 

The end result is not as important as the process of getting it. Though the team did 
not expect to win and entered the competition with an altitude to experience and 
understand rocket design and rocketry better, the results was nonetheless a good one for 
first-timers like us. There were many lessons learnt from the success and failure of some 
teams. By analyzing and discussing the causes of failure with the other teams, it provided 
much insight into the considerations that went into the design and some of the cardinal 
errors that were made. Using the example of a rocket that ‘exploded’ mid-air, their over- 
ambitious fin design, which reduced too much of the control surface for normal flight and 
stability, it provided a good lesson for our team as well in our future design. In future 
participation, more prediction tools and program should be made available for team 
members to analysis the aerodynamics and flight performance. The project may be part of 
the NPS missile design course, culminating in the actual launch of the rocket; all these 
within 2 quarters, It can help to reinforce some of the fundamental principles of 
aerodynamics and rocket design. In summary, this was a worthwhile project that 
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combined learning, application of knowledge and great fun into one, and certainly helped 
in providing a better grasp of fundamentals in missile design. 
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APPENDIX C: INDUCED DRAG PROGRAMME DEVELOPED 


1. ZLDragC_mod.m 


function [CD] = ZLDragC_mod(M,load_f,mass, rho,vel, Sref) 

% Written by LTC Weng Wai Leong, Naval Postgraduate School, Dec 2009 

% This function interpolates the known drag polars for various Mach 
number to determine drag coefficient. Input variables are Mach number, 

% load factor, mass, density, velocity, referenece area the interceptor 
or ballistic missile (BM) to be utilized 

% in the BMFlight3.m, SMTrajectory and SMFlight3.m programs (written by 
LT Lukacs) for the calculation of forces acting on the 
interceptor/Ballistic Missile. 

% The CL and Cd data are input in Excel files. Data obtained from Prof 
Yakimenko 

% Variable List 

% M = mach number of interceptor/BM 
% load_f= load factor in the normal plane 
% mass = mass of interceptor/BM 
% rho = density of air 
% vel = velocity of interceptor/BM 
% Sref = reference area of interceptor/BM 

global CL_data Cd_data 

Cd_curve = interpl(Cd_data(:,1),Cd_data(:,2:11),M, 'nearest' , 'extrap' ); 
%interpolate the Cd data in the array for the new Mach number 

CL_curve = interpl(CL_data(:,1),CL_data(:,2:11),M, ! nearest' , ’extrap 1 ); 
%interpolate the CL data in the array for the new Mach number 

p = polyfit(CL_curve,Cd_curve,2); % generate a new drag polar (CL_Cd 

curve) for new Mach number using 2nd order polynominal function 

k=p (1); 

CL0=0; 

CD0=p(3)-k*CL0 A 2; 

CL=(2*load_f*mass*9.81)/(rho*vel A 2*Sref); 

if abs(CL)>5, CL=5*sign(CL); end 

CD = CDO + k*(CL-CLO) A 2; 
return 
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APPENDIX D: MODELING PROGRAM CODE 


A. 3DOF INTERCEPTOR 
1. SMFlight3.m 


%% This is the Main M-File for the simulation. Run this file with all 
the required function files in the same current directory. 

% Written by LT John A. Lukacs IV, Naval Postgraduate School, June 2006 
% Corrected by O.Yakimenko, August 2007, November 2009 

% This script develops and tracks the flight path of the interceptor 
% missile. For the first ten seconds it integrates a series of 
% acceleration commands to simulate a vertical launch and tip over. 
Upon 

% activation of the guidance law, it sends the known values to the 
guidance 

% law and receives back the future time history of the optimal flight 
path. 

% This script then implements that optimal path. It updates the final 
% conditions and recalculates the optimal flight path at an interval of 
10 

% seconds. The script calls BRFlight3, SMParams3.m, SMDrag.m,STatmos.m, 

% and SMGuidance.m 


%% List o 
% Acc_SM 
% A11SM 
% alt 
% CD 
% count 
interval 
% dist 
% Drag 
% Forces, 
% g 


f variables 

= index-based vector of acceleration values 
= index based vector of all interceptor values 
= altitude 
= drag coefficient 

= counting variable to determine guidance law update 


SM 


= cumulative distance travelled 
= total drag force 

= index-based vector of force values 
= gravitational force, based 


on WGS-84 value of 


gravitational 

o, 

% m_i 

% Model_SM 
% MV 


attraction and altitude 
interceptor mass 

index-based vector of internal values 

= interceptor speed in Mach (relative to local speed of 


sound) 

% N1,N2 
length 
% num_SM 
% nx,ny,nz 
respectively 
% nx(y,z)_Op 
% path 
% Pos_SM 
% press 


= variables used to ensure optimal vectors are the same 

number of interations conducted (used for plotting) 

= axial acceleration command, body frame x, y and z, 

index based vector of the optimal flight path values 
returned time history of the optimal path 
index-based vector of position values 
local atmospheric pressure 
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% psi 
% psidot 
% psidot_Op 
% px,py,pz 
% px(y,z)_Op 
% py_op 
% pz_Op 

% q 
% Re 
% ro 
% Sref 
% state 
% t 

% target 
% temp 
% tgo 
% th 
% th_Op 
% thdot 
% thdot_Op 
% Thrust 
% time_Op 
% time_SM 
% update 
% V 

% V_Op 
% Vdot 
% Vdot_Op 
% Vel_SM 
% w 

interceptor 


heading angle 

rate of change of heading angle 

index based vector of the optimal flight path values 

x, y and z components of position 

index based vector of the optimal flight path values 

index based vector of the optimal flight path values 

index based vector of the optimal flight path values 

counting variable (used in another program) 

WGS-84 Earth's radius 
local atmospheric density 

planar reference area (for drag calculations) 
the state of the interceptor 
global time 

the state of the rocket 
local atmospheric temperature 
time to go to intercept 
flight path angle 

index based vector of the optimal flight path values 
rate of change of flight path angle 

index based vector of the optimal flight path values 
thrust generated by interceptor motor 

index based vector of the optimal flight path values 
index-based vector of time values 
number of updates to the guidance law conducted 
velocity of the interceptor 

index based vector of the optimal flight path values 
rate of change of the velocity 

index based vector of the optimal flight path values 
index-based vector of velocity values 

= number of 0.5s steps between the lunches of traget and 


close all, clear all, clc 


global Re alphag path % shared by SMFlight3, SMGuidance & 

SMTrajectory 

global q states % shared by SMFlight3 & SMGuidance 

global Pos_BR Pos_SM Npol Npt kpolar weight90 % ... by SMFlight3 & 

SMTrajectory 

%% Computing the flight path of a Ballistic Missile in a gravity turn 
BRFlight3 

%% Initializing variables for an Interceptor 
t=0; dt=0.5; q=0; w=120; % 60 sec detection time 
Nupd=30; count=Nupd; update=0; 

Npt=100; % the number of points the optimal trajectory is 

computed at 

Npol=8; % number of coefficients in approximating 

polynomials 

kpolar=0; 

prompt = { 'Northing wrt to BM launch point, km' , ... 

'Westing wrt to BM launch point, km',... 
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'Weighting coefficient for impact angle'}; 
dlg_title = 'Enter Interceptor launch coordinates'; 
num_lines = 1; 
def = { '80', '60', '1000' }; 

answer = inputdlg(prompt,dlg_title,num_lines, def ,' on ') ; 
px =str2num(answer{1}) *1000; 

py =str2num(answer{2}) *1000; 

weight90=str2num(answer{3}) ; 

% px=-100000/20; 

% py=100000/2; 

% weight90=100; 
pz=Re; 

px_old=px; py_old=py; pz_old=pz; 
dist=0; 

psi=atan2(Pos_BR(12 0, 2 ) -py,Pos_BR(120,1)-px) ; psi_old=psi; 
th=90*pi/180; th_old=th; 

V=l; V_o1d=V; 

%% Computing Interceptor flight path 
for i=l: 20%1000 
t=t+dt; 

%% Boost phase (vertical launch), t<10s 
if t<10 

% Speed, Mach number 
alt=norm([px;py;pz])-Re; 

[ro,press,temp]=STatmos(alt); 

MV = V/sqrt(1.402*287.053*temp); 

% Forces 

g=3.986004418e14/norm([px;py;pz]) A 2; 

[m_i,Sref] = SMParams3(t); 

CDtable = SMDrag(MV); 

if t<=6. Thrust = 206000; psidot=0; 
else Thrust = 95300; psidot=0; 

end 

ny = V/g*cos(th)^psidot; 
nz = V/g*thdot + cos(th) ; 
nn = sqrt(ny A 2+nz A 2) ; 

CL=(2*nn*m_i*g)/(ro*V A 2^Sref) ; 

CD = CDtable(1) + kpolar*CL A 2; 

Drag = ro*V A 2^CD^Sref/2 ; 
nx=(Thrust-Drag)/m_i/g; 

alphag=180/pi^nn^m_i^g/(ro*V A 2*Sref/2)/13; 

% Kinematics 
Vdot=g*(nx-sin(th)) ; 
psidot=ny^g/V/cos(th); 
thdot=g^(nz-cos(th))/V; 

% Collecting variables 
time_SM(i,1)=t; 


thdot=0; 

thdot=- 


0.075; 
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Forces_SM(i,1)=nx; 
Model_SM(i,1)=V; 
Model_SM(i,3)=th; 
Model_SM(i,5)=psi; 
Pos_SM 
Pos_SM 
Vel_SM 
Vel_SM 
Vel SM 


1.1) =px; 
i, 4)=dist; 

1.1) =V*cos(th)*cos(psi) ; 
i, 2)=V*cos(th)*sin(psi) ; 
i,3)=V*sin(th); 


Forces_SM(i, 2 ) =ny; 
Model_SM(i, 2)=Vdot; 
Model_SM(1,4)=thdot; 
Model_SM(i, 6)=psidot; 
Pos_SM(i, 2)=py; 


Acc_SM(i, 1)=Vdot*cos(th)*cos(psi)-V*cos(th )* 
-V*sin(th)*cos(psi)*thdot; 

Acc_SM(i,2)=Vdot*cos(th)*sin(psi)+V*cos(th )* 
+V*sin(th)*sin(psi)*thdot; 

Acc_SM(i,3)=Vdot*sin(th)+V*cos(th)*thdot; 


Forces_SM(i,3)=nz; 


Pos_SM(i,3)=pz; 


sin(psi)^psidot ... 
cos (psi)*psidot .. . 


% Euler integration 
V=V_o1d+Vdot *dt; 
psi=psi_old+psidot*dt; 
th=th_old+thdot*dt; 
px=px_old+V*cos(th)*cos(psi)^dt; 
py=py_old+V*cos(th)^sin(psi)*dt; 
pz=pz_old+V*sin(th)*dt; 


dist=(dist+abs(norm([px-px_old;py-py_old;pz-pz_old]))); 

V_old=V; psi_old=psi; th_old=th; 

px_old=px; py_old=py; pz_old=pz; 

else 

%% Optimal guidance, t>10s 

if count==Nupd & update==0 % Recomputing the trajectory every 

Nupd cycles 

update=update+l; 

fprintf( 1 Starting Interceptor’’s Guidance Update 

#%2.Of\n’ ,update) 

state=[px;py;pz;V;th;psi;Vdot;thdot;psidot] ; 

target=[Pos_BR(i+w,1);Pos_BR(i+w,2);Pos_BR(i+w,3); 

Vel_BR(i+w,1);Vel_BR(i+w,2);Vel_BR(i+w,3); 

Acc_BR(i+w,1);Acc_BR(i+w,2);Acc_BR(i+w,3)]; 
path=SMGuidance(t,state,target,i); % Calling SMGuidance 

function 

if update==l; 

Nl=length(path(:,1)); N2=N1; 

else 

N2=length(path ( :,1)); 

end 

% Identifying variables 

time_Op(:,update) = [path(:,1);zeros(N1-N2,1)] ; 
px_Op(:,update)=[path(:,2);zeros(N1-N2,1)]; 
py_Op(:,update) = [path(:,3);zeros(N1-N2,1)] ; 
pz_Op(:,update) = [path(:,4);zeros(N1-N2,1)]; 

V_Op(:,update)=[path(:,5);zeros(N1-N2,1)]; 
th_Op(:,update) = [path(:,6);zeros(N1-N2,1)] ; 
psi_Op(:,update)=[path(:,7);zeros(N1-N2,1)]; 

% Vdot_Op(:,update) = [path(:,8);zeros(N1-N2,1)] ; 
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thdot_Op ( : , update) = [path ( : , 9) ; zeros (N1-N2, 1 ) ] ; 
psidot_Op(:,update)=[path( : ,10);zeros(N1-N2, 1 )]; 
nx_Op(:, update) = [path( : , 8) ;zeros(N1-N2,1)] ; 
ny_Op(:, update) = [path(:,9);zeros(N1-N2, 1) ] ; 
nz_Op(:, update) = [path( :,10);zeros(N1-N2,1)]; 
count=l; 

end 

if count==101 % added by OY 2009-10-08 

count=count-l; 
elseif count==0 
count=l; 

end 


t=time_Op(count,update); 
nx=nx_Op(count,update) ; 
ny=ny_Op(count, update) ; 
nz=nz_Op(count,update) ; 

V=V_Op(count,update) ; 

Vdot=Vdot_Op(count, update) ; 
th=th_Op(count, update) ; 

thdot=thdot_Op(count, update); 
psi=psi_Op(count,update) ; 

psidot=psidot_Op(count,update); 
px=px_Op(count,update) ; 
py=py_Op(count,update); 
pz=pz_Op(count,update) ; 


Forces_SM(i,2)=ny; 
Model_SM(i,2)=Vdot; 
Model_SM(i,4)=thdot ; 
Model_SM(i,6)=psidot; 
P o s_SM(i,2)=py; 


% Collecting variables 
time_SM(i,1)=t; 

Forces_SM(i, 1)=nx; 

Model_SM(i, 1)=V; 

Model_SM(i,3)=th; 

Model_SM(i,5)=psi; 

P o s_SM(i,1)=px; 

Pos_SM(i,4)=dist ; 

Vel_SM(i,1)=V*cos(th)*cos(psi) ; 

Vel_SM(i,2)=V*cos(th)*sin(psi); 

Vel_SM(i,3)=V*sin(th) ; 

Acc_SM(i,1)=Vdot*cos(th)*cos(psi)-V*cos(th)* 
-V*sin(th)*cos(psi)*thdot; 

Acc_SM(i,2)=Vdot*cos(th)*sin(psi)+V*cos(th)* 
+V*sin(th)*sin(psi)*thdot; 

Acc_SM(i,3)=Vdot*sin(th)+V*cos(th)*thdot ; 
Update(i,1)=update; 


Forces_SM(i,3)=nz; 


Pos_SM(i,3)=pz; 


sin(psi)*psidot ... 
cos (psi)*psidot .. . 


end 


% Time step 
count=count+l; 

dist=(dist+abs(norm([px-px_old;py-py_old;pz-pz_old]))); 

V_old=V; 

psi_old=psi; 

th_old=th; 

px_old=px; 

PY_old=py; 
pz_old=pz; 

% the end of the "if" loop 
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end % the end of the "for" loop 

A11SM= [time_SM Forces_SM Model_SM Pos_SM Vel_SM Acc_SM]; % update]; 


2. SMParam3.m 


function [SM_mass, 

Sref]=SMParams3(t) 

g, 

o 

Written by LT 

John A. Lukacs IV, Naval Postgraduate School, June 2006 

g . 

This function 

calculates the J Matrix and Mass of the intercepter. 

o, 

o 

assuming a cruciform rocket in two stages to intercept. This 

g. 

function also 

returns the reference (base) diameter of the missile. 

o 

o 

% Notes: 



o. 

The first stage 

last 6 seconds, the second stage lasts an additional 

o. 

10 seconds. 

Stage 1 (booster) seperates upon completion. Stage 2 

o 

o 

does not separate after completion 

O. ; 

% List of variables 

o 

o 

dia 

= 

reference diameter, base diameter 

o 

o 

1 

= 

length, varies by component 

o. 

p_SM_stl_fuel 

= 

density of stage 1 rocket fuel 

g. 

p_SM_st2_fuel 

= 

density of stage 2 rocket fuel 

g. 

p_SMstr 

= 

density of structural material 

g. 
o 

r 

= 

radius, varies by component 

o 

o 

ro 

= 

outer radius, varies by component 

g. 

ri 

= 

inner radius, varies by component 

g. 

SM_mass 

= 

total rocket mass 

g. 

SM_nose 

= 

total mass of nosecone section 

o 

o 

SM_stl_fcr 

= 

consumption rate of stage 1 fuel 

o 

o 

SM_stl_fuel 

= 

remaining stage 1 fuel based on time and 

g. 



consumption rate 

g. 

SM_stl_str 

= 

total mass of stage 1 structural material 

g. 

SM_stl_tfm 

= 

total mass of stage 1 fuel 

g. 

SM_st2_fcr 

= 

consumption rate of stage 2 fuel 

o 

o 

SM_st2_fuel 

= 

remaining stage 2 fuel based on time and 

g. 



consumption rate 

g. 

SM_st2_str 

= 

total mass of stage 2 structural material 

g. 

SM_st2_tfm 

= 

total mass of stage 2 fuel 

g. 

t 

= 

time 

o 

o 

th 

= 

structural thickness 

g. 

V_bo dy 


= volume of body structural material 

g. 

V_nose_str 

= 

volume of nosecone structural material 

g. 

V_nose_str0 

= 

volume of nosecone structural material. 

g. 



intermediate value 

o 

o 

V_nose_str1 

= 

volume of nosecone structural material. 

o 

o 



intermediate value 

g. 

V_stl_fuel 

= 

volume of stage 1 fuel 

g. 

V_stl_str 

= 

volume of stage 1 structural material 
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% V_st2_fuel = volume of stage 2 fuel 

% V_st2_str = volume of stage 2 structural material 

%% Structural components 
p_SMstr = 4225; 
th = .0208; 

% Mass of nose cone 
1 = .8255; 
r = 0.34/2; 

V_nose_str0 = pi* (1* ( (r A 2 + l A 2 ) / (2*r) ) A 2-l A 3/3- ( ( (r A 
(2*r))-r)*((r A 2 + l A 2)/(2*r)) A 2*asin(1/ ( (r A 2 + l A 
1 = .8255-th; 
r = 0.34/2-th; 

V_nose_strl = pi* (1* ( (r A 2 + l A 2)/(2*r)) A 2-l A 3/3-(((r A 
(2 * r))-r)*((r A 2 + l A 2)/(2 *r)) A 2*asin(1/((r A 2 + l A 
V_nose_str = V_nose_strO-V_nose_strl+pi*r A 2*th; 
SM_nose = 1.3*V_nose_str*p_SMstr; 

% Mass of Body/Warhead Section 
1 = .849; 
ro = 0.34/2; 
ri = 0.34/2-th; 

V_body = l*pi*(ro A 2-ri A 2) ; 

SM_body = V_body*p_SMstr+l15; 

% Mass of Stage 1 (Mk72 Booster) 

1 = 1.72; 
ro = 0.53/2; 
ri = 0.53/2-th; 

V_stl_str = l*pi* (ro A 2-ri A 2)+2*pi*ro A 2*th; 

SM_stl_str = V_stl_str*p_SMstr; 

% Mass of Stage 2 (Mkl04 Engine) 

1 = 2.88-2*th; 
ro = 0.34/2; 
ri = 0.34/2-th; 

V_st2_str = l*pi* (ro A 2-ri A 2)+2*pi*ro A 2*th; 
SM_st2_str = V_st2_str*p_SMstr; 

%% Fuel Components 

% Stage 1 Solid Fuel is HTPB/AP/A1 

1 = 1.72; 

ri = 0.53/2-th; 

V_stl_fuel = 0.80*l*pi*ri A 2; 
p_SM_stl_fuel = 18 60; 

SM_stl_tfm = 4 68; 

SM_stl_fcr = 468/6; 

% Stage 2 Solid Fuel is TP-H1205/6 
1 = 2 . 88 ; 
ri = 0.34/2-th; 

V_st2_fuel = 0.60*l*pi*ri A 2; 

SM_st2_tfm : 360; 

p_SM_st2_fuel = SM_st2_tfm/V_st2_fuel; 


2+1 A 2)/ ... 

2)/ (2*r)))); 


2+1 A 2)/ ... 

2)/(2*r) ) ) ) ; 
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SM_st2_fcr = 360/15; 
if t<6 

%% Stage 1 - Stage 1 Fuel is consumed. Stage 2 Fuel is not used. 
SM_stl_fuel = SM_stl_tfm - SM_stl_fcr * t; 

SM_mass = SM_nose+SM_body+SM_st2_str+SM_st2_tfm+SM_stl_str+... 

SM_stl_fuel; 
dia = 0.53; 

elseif t<21 

%% Stage 2 - Stage 1 has seperated, Stage 2 Fuel is consumed. 
SM_st2_fuel = SM_st2_tfm - SM_st2_fcr * (t-6) ; 

SM_mass = SM_nose+SM_body+SM_st2_str+SM_st2_fuel; 
dia = 0.34; 

else 

%% Stage 3 - The unpowered nosecone and Stage 2 remains. 

SM_mass = SM_nose+SM_body+SM_st2_str; 
dia = 0.34; 

end 

Sref=pi*dia A 2/4; 
return 


3. SMDrag.m 


Written by LT John A. Lukacs IV, Naval Postgraduate School, June 2006 
and renamed by Prof Oleg Yakimenko in November 2009 

% This function interpolates to determine drag coefficient based on the 
% Mach number and the boost or glide phase of the rocket to be utilized 
% in the BMFlight.m and SMFlight.m programs for the calculation of 
% forces acting on the rocket/missile. 

% The tables used were point-plotted from Prof Hutchins' ME4703 
% "Missile Flight Analysis" Class Notes 

%% Setting a List of Variables 

% BDrag = boost phase drag interpolation table 
% GDrag = glide phase drag interpolation table 
% M = mach number 

% Mach = mach number interpolation table 
%% Tables: 

Mach = [0 0.90 1.1 1.2 1.5 2.0 2.5 3.0... 

3.5 5.0 6.0]; 

BDrag = [0.1444 0.1444 0.2778 0.2778 0.2308 0.1778 0.1481 0.1296... 
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0.1185 0.1000 0.0950]; 

GDrag = [0.2461 0.2461 0.4615 0.4615 0.3615 0.2846 0.2500 0.2192... 

0.2000 0.1500 0.1300]; 

%plot(Mach,BDrag,'o—',Mach,GDrag,’+-.') 

%% Calculation of Drag Coefficient Values: 
if M > 6 
M = 6; 

end 

BDrag = interpl(Mach,BDrag,M, 'cubic ') ; 

GDrag = interpl(Mach,GDrag,M, 'cubic ') ; 

Drag=[BDrag;GDrag] ; 
return 


B. 3DOF TARGET 


1. BRFlight3.m 


o 

o 

% Complementary M-File 




o 

o 

Written by 

LT John A. Lukacs TV, Naval Postgraduate 

School, June 2006 

o 

o 

This script integrates the position. 

velocity. 

and acceleration 

values 





o. 

at each time step to determine the flight 

path of 

a 

ballistic missile 

o 

o 

in a gravity turn. The script calls BRParams3.m, 

ZLDragC.m, and 

o. 

STatmos.m 





o. 

% Setting a 

List of Variables 




o. 

acc = 

total acceleration 




o. 

alt 

altitude 




o. 

ax = 

x component of acceleration 




o. 

ay 

y component of acceleration 




o 

o 

az = 

z component of acceleration 




o 

o 

CD 

drag coefficient 




o. 

Drag 

total drag force 




o. 

dt = 

time step interval 




o. 

g 

= gravitational force. 

based 

on 

WGS-84 value of 

gravitational 





o. 


attraction and altitude 




o. 

gm = 

initial launch angle 




o. 

i = 

interval count 




o. 

m_r 

rocket mass 




o 

o 

Mspd = 

rocket speed in Mach (relative 

to local 

speed of sound) 

o. 

num_BR = 

number of interations conducted (used 

for plotting) 

o. 

nx = 

axial force 




o. 

press = 

local atmospheric pressure 




o. 

px 

x component of position 




o 

o 

PY 

y component of position 




o 

o 

pz 

z component of position 
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% Re 
% ro 
% spd 
% Sref 
% t 

% temp 
% Thrust 
% vx 
% vy 
% vz 

% Acc_BR 
% Forces_BR 
% Pos_BR 
% time_BR 
% Vel_BR 
% All_BR 


WGS-84 value for Earth's radius 
local atmospheric density 
rocket speed in m/s 

planar reference area (for drag calculations) 
time 

local atmospheric temperature 
thrust generated by motor 
x component of velocity 
y component of velocity 
z component of velocity 

index-based vector of acceleration values 
index-based vector of force values 
index-based vector of position values 
index-based vector of time values 
index-based vector of velocity values 
index based vector of all rocket values 


%% Initializing Variables 
dt=0.5; i=0; 

Re—6.378137e6; 


%% Setting Initial 
px=0 ; 
py=0; 
pz=Re; 

gm=75*pi/180 ; 
vx=cos(gm); 
vy=0; 

vz=sin(gm); 


Conditions 


%% Computing Ballistic Flight Path 
for t=0:dt:200%22%19.5 
i=i+l; 


% Speed, Mach Number 
spd=norm([vx;vy;vz]); 
alt=norm([px;py;pz])-Re; 
if alt<86000 

[ro,press,temp]=STatmos(alt); % Calling STatmos function 

else 

[ro,press,temp]^STatmos(86000); % Calling STatmos function 

end 

Mspd=spd/sqrt(1.402*287*temp); 


% Forces 

g=3.986004418e14/norm([px;py;pz]) A 2; 

[m_r,Sref]=BRParams3(t); % Calling BRParams3 

function 

CD=BRDrag(Mspd); % Calling BRMrag function 


if t<125 

Thrust=105000^ 9.81 ; 
CD=CD (1) ; 
elseif t<240 
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Thrust=29950*9.81; 

CD=CD (1) ; 

else 

Thrust=0; 

CD=CD (2) ; 

end 

Drag=ro*spd A 2*CD*Sref/2 ; 
nx=(Thrust-Drag)/m_r/g; 

% Accelerations 

g=3.986004418el4*pz/norm([px;py;pz]) A 3; 

ax=g*nx*cos(gm); 

ay=0; 

az=g*nx*sin(gm)-g; 
acc=norm([ax;ay;az]); 

% Collect Variables 
time_BR(i ,1) =t; 

P o s_BR(i r 1) =px; 

P o s_BR(i,2)=py; 

P o s_BR(i , 3)=pz; 

Pos_BR(i,4)=norm([px;py;pz]); 

Ve1_BR(i , 1)=vx ; 

Vel_BR(i,2)=vy; 

Vel_BR(i r 3)=vz ; 

Vel_BR(i,4)=spd; 

Vel_BR(i,5)=Mspd; 

Acc_BR(i , 1)=ax; 

Acc_BR(i f 2)=ay; 

Acc_BR(i,3)=az; 

Acc_BR(i,4)=acc; 

Acc_BR(i , 5)=nx; 

Forces_BR(i , 1)=Thrust; 

Forces_BR(i , 2)=m_r; 

Forces_BR(i, 3)=Drag; 

% Time Step 

px=px+dt *vx ; 

py=py+dt^vy ; 

pz=pz+dt*vz ; 

vx=vx+dt*ax; 

vy=vy+dt^ay; 

vz=vz+dt*az ; 

num_BR=length(time_BR); 

end 

%plot(time_BR( :,1), (Pos_BR(: f 3)-Re)/1000) 

A11BR= [time_BR Forces_BR Pos_BR Vel_BR Acc_BR]; 

clear CD Drag Mspd Sref Thrust acc alt ax ay az dt 

clear g gm h i m_r nx press px py pz ro spd t temp vx vy vz 
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2. BRParams3.m 


function [BR_mass,Sref,length]=BRParams(t) 

% Written by LT John A. Lukacs IV, Naval Postgraduate School, June 2006 

% This function calculates the mass of the rocket, assuming a cruciform 
% rocket in two stages plus an unpowered nosecone stage. This function 
% also returns the reference (base) diameter of the missile. 


%% Notes: 

% The first stage last 125 seconds, the second stage lasts an 
additional 

% 110 seconds. The stages seperate upon completion. 


%% Setting a List of Variables 
% BR_mass = total rocket mass 
% BR_nose = total mass of nosecone section 
% BR_stl_bt = burntime for stage 1 

% BR_stl_fcr = consumption rate of stage 1 fuel 

% BR_stl_fuel = remaining stage 1 fuel based on time and consumption 


rate 

% BR_stl_str 
% BR_stl_tfm 
% BR_st2_bt = 
% BR_st2_fcr 
% BR_st2_fuel 
rate 

% BR_st2_str 
% BR_st2_tfm 
% dia 


= total mass of stage 1 structural material 

= total mass of stage 1 fuel 

burntime for stage 2 

= consumption rate of stage 2 fuel 

= remaining stage 2 fuel based on time and consumption 

= total mass of stage 2 structural material 

= total mass of stage 2 fuel 

= reference diameter, base diameter 


% t = time 


%% Setting Structural Components 
BR_nose = 250; 

BR_st2_str = 2288; 

BR_stl_str = 9000; 

%% Setting Fuel Components 
BR_stl_tfm = 50970; 

BR_stl_bt = 125; 

BR_stl_fcr = BR_stl_tfm/BR_stl_bt; 

BR_st2_tfm = 12912; 

BR_st2_bt = 110; 

BR_st2_fcr = BR_st2_tfm/BR_st2_bt; 
if t<125 

%% Stage 1 - Stage 1 Fuel is consumed. Stage 2 Fuel is not used 
BR_stl_fuel = BR_stl_tfm-BR_stl_fcr*t; 

BR_mass = BR_nose+BR_stl_str+BR_stl_fuel+BR_st2_str+BR_st2_tfm; 
dia = 2.2; 
length = 2+14+16; 

elseif t<240; 

%% Stage 2 - Stage 1 has seperated. Stage 2 Fuel is consumed 
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BR_st2_fuel = BR_st2_tfm-BR_st2_fcr*(t-125); 

BR_mass = BR_nose+BR_st2_str+BR_st2_fuel; 
dia = 1.3; 
length = 2 + 14; 

else 

%% Stage 3 - Stage 2 has seperated, only the unpowered nosecone remains 
BR_ma s s = BR_n ose; 
dia = 1.3; 
length = 2; 

end 

Sref=pi*dia A 2/4; 
return 


3. BRDrag.m 


function Drag = BRDrag(M) 

% Written by LT John A. Lukacs IV, Naval Postgraduate School, June 2006 

% Renamed by Oleg Yakimenko, November 2009 to be used as separate drag 
file 

% for interceptor and ballistic missile instead of the common ZLDragC.m 
file 

% used previously 

% This function interpolates to determine drag coefficient based on the 
% Mach number and the boost or glide phase of the rocket to be utilized 
% in the BMFlight.m and SMFlight.m programs for the calculation of 
% forces acting on the rocket/missile. 

% The tables used were point-plotted from Prof Hutchins' ME4703 
% "Missile Flight Analysis" Class Notes 

%% List of variables 


% BDrag = boost phase drag interpolation table 
% GDrag = glide phase drag interpolation table 
% M = mach number 

% Mach = mach number interpolation table 


%% Tables: 


Mach 

[0 

3.5 

0.90 

5.0 

1.1 

6.0] ; 

1.2 

BDrag 

[0.1444 

0.1444 

0.2778 

0.2778 


0.1185 

0.1000 

0.0950] 

r 

GDrag = 

[0.2461 

0.2461 

0.4615 

0.4615 


0.2000 

0.1500 

0.1300] 

i 


1.5 2.0 2.5 

0.2308 0.1778 0.1481 
0.3615 0.2846 0.2500 


3.0. . . 

0.1296. . . 
0.2192. . . 
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%plot(Mach,BDrag,'o—',Mach,GDrag,'+-.') 

%% Calculation of drag coefficient for the boost and glide phases: 
if M > 6 
M = 6; 

end 

BDrag = interpl(Mach,BDrag,M, 'cubic ') ; 

GDrag = interpl(Mach,GDrag,M, 'cubic ') ; 

Drag=[BDrag;GDrag] ; 
return 


C. GUIDANCE ALGORITHMS 
1. SMGuidance.m 


function path=SMGuidance(time,state,target,i) 

% Written by LT John A. Lukacs IV, Naval Postgraduate School, June 2006 
% Corrected by O.Yakimenko, August 2007, October 2009 


% This function takes in the state of the interceptor and target and 
% generates an initial guess at the final conditions (position, 

% orientation angles, range, and time to intercept) through a first- 
order 

% trajectory assumption and iterative process. It then calls the 
% fminsearch function using those initial guesses. Finally, it plots 
the 

% returned optimal flight path and associated variables. 


%% List of 
% best 

o, 

% BC 
% cost 
% costs 
iteration 
% free 

o, 

% init 
% J 
% N 

% nmax 

o, 

% path 
% psi 
% psidot 
% psif 

o, 

% psit 
% Py,Pz 
% q 


variables 

= vector of the variables in the optimal path returned from 
the fminsearch function 
= boundary conditions 

= cost function value returned from SMGuidanceCost function 
= array of the value of the cost variables at each 

= variables that fminsearch can modify, specifically 
[ tau;thf;psif] 

= vector of initial estimates 
= vector of cost function variable values 
= length of the path vector (used for plotting) 

= maximum acceleration capability of the interceptor, 
altitude dependent 

= returned time history of the optimal path 
= initial interceptor heading angle 

= initial rate of change of interceptor heading angle 
= final interceptor heading angle, calculated from final 
conditions estimate 
= target heading angle 

= penalty functions on the y and z acceleration 
= variable, counting trajectory updates during intercept 
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% qq 

= variable, counting the number of iterations during 

optimization 


% ranqe 

= estimate of distance between target and interceptor 

% state 

= state of the interceptor missile, sent from SMGuidance.m 

o, 

o 

[px;py;pz;V;th;psi;Vdot;thetadot;psidot]; 


% states 

= array of the values of all processes in SMGuidance.m 

% tarqet 

= state of the rocket, sent from SMGuidance.m at i+w. 

o, 

o 

synchronizing times 


% tau_f 

= value of the virtual arc 


% tgo 

= time to go to intercept 


% th 

= initial interceptor flight path angle 


% thdot 

= initial rate of change of interceptor flight 

path angle 

% thf 

= final interceptor flight path angle 


% tht 

= target flight path angle 


% tic..toe 

= MATLAB function to track run time 


% trys 

= vector of optimal path and derivative values 


% V 

= initial interceptor velocity 


% V_f 

= final interceptor velocity 


% Vave 

= average interceptor velocity 


% Vdot 

= initial interceptor acceleration 


% xO 

= initial inteceptor position 


% xdO 

= initial interceptor velocity 


% xddO 

= initial inteceptor acceleration 


% xdf 

= final inteceptor position 


% xdt 

= current target velocity 


% xf 

= final inteceptor acceleration 


% xmult 

= ratio value (used for plotting) 


% xt 

= current target position 


% ymult 

= ratio value (used for plotting) 


% zmult 

= ratio value (used for plotting) 


global Re 

alphag path % shared by SMFlight3, 

SMGuidance & 

SMTrajectory 


global q states % shared by SMFlight3 & SMGuidance 

global qq 

thdot psidot tgo costs trys % ... by 

SMGuidance & 

SMTrajectory 


%% Counting trajectory updates during intercept 


q=q+l; 



%% Initializing Interceptor (all states) 


xO=state (1 

: 3) ; 


V=state (4) 

I 


th=state (5) ; 


psi=state ( 

6) ; 


Vdot=state 

(7) ; 


thdot=state (8); 


psidot=state (9) ; 


xdO = [V* 

cos(th)*cos (psi) ; 


V* 

cos(th)*sin (psi) ; 


V* 

sin(th)] ; 


xddO 

= [Vdot*cos(th)*cos(psi)-V*cos(th)*sin(psi)*psidot- 

V*sin(th)* 

cos(psi)*thdot; 


Vdot*cos(th)*sin(psi)+V*cos(th)*cos(psi)*psidot- 


V*sin (th)* 

sin(psi)*thdot; 
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Vdot*sin(th)+V*cos(th)*thdot] ; 

%% Initializing BM (up to the second-order derivatives) 
xt =target(1:3); 
xdt =target(4:6); 
xddt=target(7:9) ; 

%% Estimating time-to-go 
tgol=100; delta=5; 
vmaxhyp=2 60 0; 
while delta>l 
if tgol>20 

V_f=vmaxhyp-10*(tgol-20); 

Vave=(20*vmaxhyp/2+(tgol-20)*(vmaxhyp+V_f) / 2 ) /tgol; 

else 

V_f=tgol*vmaxhyp/20; 

Vave=V_f/2; 

end 

xf : =xt+xdt*tgol; 

tgo2=sqrt((xf(l)-xO(1)) A 2+(xf(2)-xO(2)) A 2+(xf(3)-xO(3)) A 2) ... 

/(norm(xdt)+Vave) ; 
delta=abs(tgo2-tgol) ; 
tgol=(tgol+tgo2)/2; 

end 

tgo=tgol; 

%% Initializing optimization 

range=sqrt((xf(1)-xO(1)) A 2+((xf(2)-xO(2))) A 2+((xf(3)-xO(3))) A 2); 

%fprintf('Trajectory update # %2.0f \n',q) 

fprintf( '\nSlant range to target: %5.1fkm \n’ , range/10 A 3) 

tau_f =0.00045^range-1000*(q-1); % guess on tau_f 

fprintf(' Guess on the virtual arc length: %5.1f \n' ,tau_f) 
fprintf(' Guess on the time-to-go: %5.1fkm \n',tgo) 

predicted_xdt=xddt^tgo; 

tht =atan2(predicted_xdt ( 3 ), norm(predicted_xdt(1:2))); 
psit=atan2(predicted_xdt(2) , predicted_xdt(1)) ; 
thf =0;%-tht; % guess on thf 

psif =psi; %psit+pi; % guess on psif 

free =[tau_f;thf;psif;.1;1;1;-0.001;-0.001]; 

BC=[xO;xdO;xddO;time;xt;xdt; xddt] ; 

%% Searching for the minimum performance index 
qq=0; % counting iterations to converge 

tic 

options=optimset( 'Maxiter ',100, 'Tolfun',1, 1 ToIX' ,1); 

best = fminsearch(@(x) SMTrajectory(x,BC),free,options); % 

Optimization 

tcpu=toc; 

fprintf (' \nlt took %6.0f interations to converge\n' ,qq) 

fprintf (' Elapsed time is %6.1f seconds\n' ,tcpu), 

fprintf(' Combined performance index is %5.If\n ', costs(end,1)) 

fprintf( ' including: tau_f=%5.1f, t2go=%5.1fs, ImpAngle=%4.If 

off,\n' , . . . 

costs(end,2),costs(end,3),180/pi*costs(end,4)) 
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fprintf( ' Pny=%5.1f and 

Pnz=%5.If\n' , costs(end,5), costs(end, 6)) 

%fprintf(' Dtgo=%5.1f and 

Altfine=%5.If\n',costs(end,7),costs(end,8)) 

tau_f=best(1) ; 
thf =best(2); 
psif =best(3); 

[best (4);best(5);best(6);best(7);best(8)] ; 

V_f=path(end, 5) ; 
xf=path(end,2:4); 

fprintf([ 1 Impact occurs at Altitude of %4.1fkm, Northing=%4.Ifkm, and 

i 

’Westing=%4.Ifkm\n' ] ,[xf(3)-Re xf(1) 

xf (2) ]/10 A 3) 

%fprintf('with interceptor’’s speed of %4.1f km/s\n',V_f/10 A 3) 

xdf=[V_f*cos(thf)*cos(psif) ; 

V_f*cos(thf)*sin(psif); 

V_f*sin(thf) ] ; 

%% Plotting results 

%{ 

xmult=20000/(norm(xdO)+norm(xdt)); 
ymult=20000/(norm(xdO)+norm(xdt)); 
zmult=20000/(norm(xdO)+norm(xdt)); 
nmax=40+(40-10)/(0-50000)*(path(:,4)-Re); 

figure('NameBird-eye view') % Bird-eye view 

plot3 (path(:,2)/10 A 3,path(:,3)/10 A 3, (path(:,4)-Re)/10 A 3, ’- 
.b',’Linewidth',2) 
hold on, grid 

plot3(10 A -3*[xO(1)-xmult*xd0(1);xO(1)+xmult*xd0(1)],... 

10 A -3*[xO(2)-ymult*xd0(2);x0(2)+ymult*xd0(2)],... 

10 A -3*[xO(3)-Re-zmult*xd0(3);x0(3)- 
Re+zmult*xd0(3)], 1 c 1 , 1 Linewidth' , 2) 

plot3(xf(1)/10 A 3,xf(2)/10 A 3, (xf(3)-Re)/10 A 3, ’pr', 'Linewidth',2) 
plot3(10 A -3*[xf(1)-xmult*xdt(1);xf(1)+xmult*xdt(1)],... 

10 A -3*[xf(2)-ymult*xdt(2);xf(2)+ymult*xdt(2)],... 

10 A -3* [xf (3) -Re-zmult*xdt (3);xf (3)- 
Re+zmult*xdt(3) ] , ' r ' , 'Linewidth' , 2) 

plot3(xO(1)/10 A 3,xO(2)/10 A 3,(xO(3)-Re)/10 A 3, 1 *b 1 ,'Linewidth',5) 
plot3((xO(1)+xmult*xd0(1))/10 A 3,(xO(2)+xmult*xd0(2))/10 A 3,... 

(xO (3) - 

Re+xmult*xd0(3))/10 A 3,’ A c’,'linewidth 1 ,2) 

plot3((xf(1)+xmult*xdt(1))/10 A 3,(xf(2)+xmult*xdt(2))/10 A 3,... 

(xf (3)- 

Re+xmult*xdt(3))/10 A 3, ,A r’, 'linewidth' , 2) 

hl=legend('Intercept trajectoryInterceptor''s velocity vector’,... 

’Impact point’,'BM’’s velocity vector at 

intercept' , ’Location’ , ’Best'); 
set (hi, 'FontSize’,8); 

xlabel('Northing (km)'), ylabel('Westing (km)’), zlabel('Altitude 
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(km)') 

%title ('Interception Geometery', 'Fontsize',10) 
view(-102,8) 
axis equal 

figure('Name' , 'Combined PI and Virtual arc') % Performance index & 

Virtual arc 
subplot(211) 

semilogy(costs(:,1)/costs (1,1), 'h-.'), grid 

xlabel('Iteration'), ylabel('Relative PI (PI_i/PI_l') 

xlim([l qq]) , axis 'auto y' % ylim([0 2]) 

subplot(212) 

plot (costs(:,2),'h-.'), grid 

xlabel('Iteration'), ylabel('Length of virtual arc, \it\tau_f') 
xlim([1 qq]) 

figure('Name','Delta Time-to-go and Altitude violation') % Dtgo & 
Altitude fine 
subplot (211) 

plot(costs(:, 7),'h-.') , grid 
xlim([1 qq]) 

xlabel('Iteration') , ylabel('\Delta \itt_{go} \rm(s)') 
subplot(212) 

plot(costs(:,8)/10 A 3, 'h-. ') , grid 

xlabel('Iteration'), ylabel('Alt. violation, (km)') 
xlim([1 qq]) 

figure('Name' , 'Impact angle and Time-to-go') % Impact angle & Time- 

to-go 

subplot(211) 

plot(real(180/pi*acos(costs (:,4))),'h-.'), grid 
axis([1 qq 60 90] ) 

xlabel ('Iteration'), ylabel('Impact angle ( A o)') 
subplot(212) 

plot (costs ( :,3), 'h-. '), grid 

xlabel('Iteration') , ylabel('Time-to-go, \itt_{go} \rm(s) ') 
xlim([1 qq]) 

figure('Name','G-load factors') % G-load constraints 

subplot(211) 

plot (path(:,1),path(:,9), '-b. ') , grid 
hold on 

plot (path(:,1),path(:,10) , '—g. ' ) 
plot(path(:,1),nmax(:),'r','Linewidth',2) 
plot(path(:,1),-nmax(:),'r','Linewidth',2) 
hl=legend('n_y','n_z','Dynamic constraints',2); 
set(hi,'FontSize',8); 

xlabel('Time (s)'), ylabel('Load factor (g)') 
subplot(212) 

plot(costs(:,5)/10 A 7,'-b.','Linewidth',2), grid 
hold on 

plot(costs(:,6)/10 A 7,'—g.','Linewidth',2) 
xlabel('Iteration') , ylabel('Relative penalty') 
hl=legend('n_y penalty','n_z penalty','Location','Best'); 
set(hi,'Fontsize',8); 
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xlim([1 qq]) 

figure('Name', ’Lambda and Tau profile') % Lambda and tau 

subplot(211) 

plot (path(:,1),path(:,end), ' . ') , grid 
xlabel('Time (s)'), ylabel('\it\lambda') 
subplot(212) 

plot (path ( : , 1) ,path(:,end-1),'.'), grid 
xlabel('Time (s)'), ylabel( 1 \it\tau') 

figure('Name', 'Speed and Angle of attack profile') % SM speed and Angle 
of attack 
subplot (211) 

plot(path( : , 1) ,path(:,5),'.') ; grid 
xlabel('Time (s)'), ylabel('Speed, V (m/s)') 
subplot(212) 

plot (path(:,1),alphag, ' . ') , grid 

xlabel('Time (s)'), ylabel('Angle of attack ( A o)') 

figure('Name', 'Euler angles profile') % SM Euler angles 

subplot(211) 

plot (path( : , 1) ,path(:,6)*180/pi, ' . ' , 'Linewidth',2), grid 
hold on 

plot(path(end,1),thf*180/pi,'ro') 
ylim([-30 90]) 

xlabel('Time (s)'), ylabel('\theta ( A o)') 
subplot(212) 

plot (path {:,!), path(:,7)*180/pi, 'g. ', 'Linewidth',2), grid 
hold on 

plot(path(end,l),psif*180/pi,'ro') 
ylim([-180 180] ) 

xlabel('Time (s)'), ylabel('\psi ( A o)') 

%} 

%% Creating results structure 

states{q,1}=path; 

states{q,2}=BC; 

states{q,3}=free; 

states{q,4}=best; 

states{q,5}=costs; 

return 


2. SMGuidanceCost.m 


function [cost,J,Py,Pz]=SMGuidanceCost(free,const) 

% Written by LT John A. Lukacs IV, Naval Postgraduate School, June 2006 

% This function calculates the cost of the proposed trajectory returned 
% from the SMTrajectory.m function based on the optimization parameters 
% and penalty parameters defined herein. This is a sub-funtion of the 
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% SMGuidance.m function's fminsearch. This cost value is used to 
% determine whether the proposed trajectory is optimal. The trajectory 
% that returns the minimum value of J is the optimal function. 


%% Variable 
% calccost 

o, 

o 

% const 

o, 

% cost 
function 
% costs 
% dist 
% free 

o, 

% init 
% J 
% N 

% nmax 

o, 

% nx 
% ny 
% nz 
% path 

o, 

nz ' ] 

% psi 
% psidot 
% psif 

o, 

% Py 
% Pz 

% qq 
% t 

% tau_f 
% tgo; 

% th 
% thdot 
% thf 
% time 
% V 
% V_f 
% Vdot 
% X 
% xO 
% xdO 
% xddO 
% xdf 
% xdt 
% xt 


List 

= a global variable of the value of the J function (used 
for plotting) 

= variables that fminsearch cannot modify, including system 
contraints, specifically [xO;xdO;xddO;time;xt;xdt;init] 

= cost function value returned from SMGuidanceCost.m 

= vector of values of the cost variables at each iteration 
= cumulative distance travelled 

= variables that fminsearch can modify, specifically 
[tau;tgo;thf;psif] 

= vector of initial estimates 
= vector of cost function variable values 
= length of the path vector (used for plotting) 

= maximum acceleration capability of the interceptor, 
altitude dependent 

= axial acceleration command, body frame x 
= axial acceleration command, body frame y 
= axial acceleration command, body frame z 

= returned time history of the optimal path, specifically 
[time' X(l:3,:)' V' th' psi' Vdot' thdot' psidot' nx' ny' 

= initial interceptor heading angle 

= initial rate of change of interceptor heading angle 
= final interceptor heading angle, calculated from final 
conditions estimate 

= penalty function on the y acceleration 

= penalty function on the z acceleration 

= counting variable 

= current time 

= value of the virtual arc 

= time to go to intercept 

= initial interceptor flight path angle 

= initial rate of change of interceptor flight path angle 
= final interceptor flight path angle 
= optimal path time history 
= initial interceptor velocity 
= final interceptor velocity 
= initial interceptor acceleration 

= the optimal path time history in cartesian coordinates 
= initial inteceptor position 
= initial interceptor velocity 
= initial inteceptor acceleration 
= final inteceptor velocity 
= current target velocity 
= current target position 


[path]=SMTrajectory(free, const) ; 


global calccost costs q qq tgo trys 
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% Initialize Variables 

qq=qq+l; 

dist=0; 

Re=6.378137e6; 

%% Identify Variables 
time=path(: , 1 ); 

X=path(:,2:4); 

V=path (:, 5); 
th=path (:, 6); 
psi=path {:,!)} 

Vdot=path(: , 8) ; 
thdot=path(: , 9) ; 
psidot=path( :,10); 
nx=path(: , 11) ; 
ny=path {:, 12 ) ; 
nz=path(:,13) ; 

N=length(path {:,!)); 

tau_f=free(1) ; 
tgo=free(2) ; 
thf=free(3); 
psif=free(4); 

x0=const (1:3); 
xd0=const (4:6); 
xdd0=const (7:9); 
t=const (10); 
xt=const (11:13); 
xdt=const(14:16); 
init=const (17:19); 

V_f=path(N,5); 

xdf=[V_f*cos(thf)*cos(psif); 

V_f*cos(thf)*sin(psif); 

V_f*sin(thf) ] ; 
tgo=path(N,1) ; 

for i=2:1:N 

dist=dist+abs(norm([X(i,1)-X(i-1,1);X(i f 2)-X(i-l f 2);X(i f 3)-X(i- 
1,3);])); 
end 

nmax=40+(40-10)/(0-50000)*(path(:,4)-Re); 

%% Calulate Cost of the chosen trajectory 
J=[ tau_f; 
tgo; 

100*abs(dot(xdf , xdt)/norm(xdf)/norm(xdt))]; 

Py=sum(max(0 , abs(ny)-nmax) . A 2 ) ; 

Pz=sum(max(0 , abs(nz)-nmax). A 2); 

cost=0.33*ones(1^3)*J+norm([Py;Pz]); 
costs(qq,1:6) = [cost;J;Py; Pz; ] ; 
calccost=cost; 
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return 


3. SMTrajectory.m 


function cost=SMTrajectory(free,BC) 

% This function computes a candidate trajectory and associated cost 
based on 

% the vector of varied parameters "free” and boundary conditions "BC". 

% This is a sub-funtion of the SMGuidance.m function's fminsearch. 

% This function creates a 7th order set of equations and evaluates that 
set at 

% the boundary conditions supplied by the inputs. It then calculates 
the time 

% history of all the flight vehicle variables, including controls and 
reactions, 

% necessary to develop that flight path. A plot command set at the end 
of this 

% function will plot a chart of the iterations at the end of run if 
desired. 

% Finally, this function calculates the cost of the candidate 
trajectory 

% combining the value of the performance index and penalties. This cost 
value is 

% used to determine whether the proposed trajectory is optimal. (The 
trajectory 

% that returns the minimum value of the cost is the sub-optimal one.) 

% O.Yakimenko, Naval Postgraduate School, November 2009 


%% List of variables 

% A,Ax,Axp,Axpp,Axppp = cell matrices of coefficients of a candidate 
reference 

% trajectory and their derivatives wrt virtual 


arc 

% BC = boundary conditions, specifically 

[xO;xdO;xddO;time;xt;xdt;xddt] 

% Cx,Cxp,Cxpp,Cxppp = coefficients of a candidate reference 

trajectory and 

their derivatives wrt virtual arc 
= tau step value 
= time step value 

= variable parameters, specifically [tau;thf;psif] 

= gravitational force 
= lambda, virtual speed 

= first-order derivative of lambda wrt to virtual arc 
= maximum acceleration capability of the interceptor, 
altitude dependent 

% nX,nXp,nXpp,nXppp = norm of reference trajectory and its 


% dtau 
% dtime 
% free 

% g 
% L 
% Lp 
% nmax 

o, 

o 
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derivatives 
% nx 
% ny 
% nz 
% path 

o, 

% psi 
% psidot 

% qq 
% t 
% tau 
% tau_f 
% tgo 
% th 
% thdot 
% thp 

virtual arc 
% time 
% trys 
plotting) 

% V 

% Vdot 
% Vp 

% X,Xp,Xpp, 


= axial acceleration command, body frame x 

= axial acceleration command, body frame y 

= axial acceleration command, body frame z 

= returned time history of the optimal path, specifically 
[time' X(l:3,:)' V 1 th ’ psi’ nx' ny 1 nz 1 tau' L 1 ] 

= interceptor heading angle 

= rate of change of interceptor heading angle 
= variable counting the number of iterations 
= global current time 
= virtual arc 

= length of the virtual arc 
= time to go to intercept 
= interceptor flight path angle 

= rate of change of interceptor flight path angle 

= first-order derivative of flight path angle wrt to 


= optimal path time history time=[0;tgo] 

= collection of norms [X nX Xp nXp Xpp nXppp] 


(used for 


= velocity 
= acceleration 

= irst-order derivative of velocity wrt to virtual arc 
Xppp = reference trajectory and its derivatives wrt 


tau 

% xO,xf = initial 

% xdO,xdf = initial 
% xddO,xddf = initial 
% xp,xpp,xppp 
virtual domain 
% xt,xdt,xddt 
time=0 


and final inteceptor position 
and final interceptor velocity 
and final inteceptor acceleration 

= boundary conditions (at 0 and f) in the 

= target position, velocity and acceleration at 


global Re alphag path % shared by SMFlight3, SMGuidance & 

SMTrajectory 

global Pos_BR Pos_SM Npol Npt kpolar weight90 % ... by SMFlight3 & 

SMTrajectory 

global qq thdot psidot tgo costs trys % ... by SMGuidance & 

SMTrajectory 

%% Counting iterations to converge 

qq=qq+l; 

tgoold=tgo; 

%% Assigning variables 
tau_f=free(1); 
thf =free(2); 
psif =free(3); 

t=BC(10); % current time in the SM frame to compute it's stage (thrust 

and drug) 

xO=BC (1:3); 
xdO=BC(4:6); 
xddO=BC(7:9); 

V(1)=norm(xdO); 
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Vdoti=norm(xddO) ; 


L (1)=1; %V(1); 
Lpi=0; %Vdoti/L(1) 


xt=BC(11:13); 
xdt=BC(14:16) ; 
xddt=BC(17:19) ; 


% Target's coordinates at the moment of detection 
% Target's velocities at the moment of detection 
% Target's accelerations at the moment of detection 


V(Npt)=2600-1/3*(tgo-20); % SM velocity estimate at impact 

Vdotf=-0.3; %-5.9578; % SM acceleartion estimate at impact 


L(Npt)=1; %V(Npt); 

Lpf=0; %Vdotf/L(Npt) ; 

xf : =xt+xdt*tgo+0.5*xddt*tgo A 2; 
xdf=[V(Npt)*cos(thf)*cos(psif); 
V(Npt)*cos(thf)*sin(psif) ; 
V(Npt)*sin(thf)]; 


thdotf=free(7); psidotf=free(8); 

xddf=[Vdotf*cos(thf)*cos(psif)-V(Npt)*cos(thf)*sin(psif)*psidotf- ... 
V(Npt)*sin(thf)*cos(psif)*thdotf; 

Vdotf*cos(thf)*sin(psif)+V(Npt)*cos(thf)*cos(psif)*psidotf-... 
V(Npt)*sin(thf)*sin(psif)*thdotf; 

Vdotf*sin(thf)+V(Npt)*cos(thf)*thdotf] ; 


%% Converting boundary conditions ito the virtual domain 
xp0=xd0/L(1) ; 

xpp0=(xddO-xdO*Lpi)/L(1) A 2; 
xppp0=[free(4);free(5);free(6)]; 

xpf=xdf/L(Npt); 

xppf=(xddf-xdf*Lpf)/L(Npt) A 2; 

xpppf~[0;0; 0] ; 

%% Calculating polynomials' coefficients and reference trajectories 
dtau =tau_f/(Npt-1); 
tau =linspace(0 , tau_f , Npt); 
for i=l:3 
A{i}=[x0(i); 

xpO (i); 
xppO (i); 
xpppO (i); 

(-16*xppp0(i)-4*xpppf(i))/tau_f+(- 
120*xpp0(i)+60*xppf(i))/tau_f A 2+... 

(-360*xpf(i)-480*xp0(i))/tau_f A 3+(840*xf(i)- 
840*x0 (i))/tau_f A 4; 

(60*xppp0(i)+30*xpppf(i))/tau_f A 2+(600*xpp0(i)- 
420*xppf(i))/tau_f A 3+... 

(2340*xpf(i)+2700*xp0(i))/tau_f A 4+(5040*x0(i)- 
5040*xf(i))/tau_f A 5; 

(-80*xppp0(i)-60*xpppf(i))/tau_f A 3+(780*xppf(i)- 
900*xpp0(i))/tau_f A 4+... 
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(-4080* xp f(i)-4320* xp 0(i))/tau_f A 5+(- 
8400*x0(i)+8400*xf(i))/tau_f A 6; 


(35*xppp0(i)+35*xpppf(i))/tau_f A 4+(420*xpp0(i)- 
420*xppf(i))/tau_f A 5+... 



(2100*xpf(i)+2100* 

xpO(i))/tau_f A 

6+ (4200 

*x0(i) 

- 


4200*xf(i))/tau_f A 7]; 







Ax{i} 

])*A{i}; 

=diag ( [ 1, 1, 

1/2, 

1/6, 

1/24, 

1/60, 

1/120, 

1/210 

Axp{i} 

])*A{i}; 

=diag ( [ 0, 1, 

1, 

1/2, 

1/6, 

1/12, 

1/20, 

1/30 

Axpp{i} 

])*A{i}; 

=diag ( [ 0, 0, 

1, 

1, 

1/2, 

1/3, 

1/4, 

1/5 

Axppp{i} 

=diag([ 0, 0, 0, 

1, 

1, 

1, 

1, 

1 ])*A{i}; 


Cx(i, :) -Ax{i}([Npol:-1:1]); 

Cxp(i,:) =Axp{i}([Npol:-1:2]); 

Cxpp(i,:) =Axpp{i}([Npol:-1:3]); 

Cxppp(i,:)=Axppp{i}([Npol:-1:4]); 

X(i,:) =polyval(Cx(i,:),tau); 

Xp(i, :) =polyval(Cxp(i,:),tau); 

Xpp (i, : ) =polyval(Cxpp(i, : ) ,tau); 

Xppp (i, :) =polyval(Cxppp(i, :),tau); 

end 

%% Computing the states 
xpl2=Xp(l,:). A 2+Xp(2,:). A 2; 
th =atan2(Xp(3,:),sqrt(xpl2)); 
psi =atan2(Xp(2, :) ,Xp(1, :)) ; 

thp =(Xpp(3,:).*xpl2- 

Xp(3,:).*(Xp(1,:).*Xpp(1,:)+Xp { 2 ,:) .*Xpp(2,:)))./.. . 

sqrt(xpl2) ./ (xpl2+Xp ( 3 , :) . A 2) ; 

psip =(Xp(1,:).*Xpp { 2 ,:) -Xpp(1,:).*Xp(2,:))./xpl2; 
time(1)=t; 

g = 3.986004418el4/norm(X(1:3,1)) A 2; 

nx(1) = Vdoti/g+sin(th(1)) ; 

ny(l) = V(1)/g*cos(th(1))*psidot; 
nz(1) = V(1)/g*thdot + cos(th(1) ) ; 

[ro,press,temp]=STatmos(norm(X(:,l))-Re); 

[m_i A Sref]=SMParams3(time(1)); 

alphag(1)=180/pi*sqrt(ny(1) A 2+nz(1) A 2)*m_i*g/(ro*V(l) A 2*Sref/2)/13; 
for j=2:Npt; 

g=3.986004418el4/norm(X(l:3 f j)) A 2; 

if norm(X(:,j))-Re<86000 , [ro,press,temp]=STatmos(norm(X(: A j))- 

Re) ; 

else [ro,press,temp]=STatmos(86000); 

end 

MV=V(j — 1)/sqrt(1.402*287.053*temp) ; 

CDtable^SMDrag(MV); 

if time(j-l)<20 Thrust=95300; CD0=CDtable(1); 

else Thrust=0; CD0=CDtable(2); 

end 

[m_i,Sref]=SMParams3(time(j — 1) ); 


95 



CL=(2*sqrt(ny(j-l) A 2+nz(j—1) A 2)*m_i*g)/(ro*V(j-1) A 2*Sref); 
CD = CDO + kpolar*CL A 2; 

Drag=ro*V(j—1) A 2*CD*Sref/2; 
nx(j)=(Thrust-Drag)/m_i/g; 

V(j)=V(j—1)+g*(nx(j —1)-sin(th(j —1)))/L(j —1)*dtau; 

ddist = sqrt( (X(1, j)-X(1,j-1) ) A 2+(X(2, j)-X(2, j-1) ) A 2+(X(3,j)-X(3,j- 
1) ) A 2) ; 

dtime=2*ddist/(V(j)+V(j-1)); 

L(j)=dt au/dtime; 

ny(j)=V(j)/g*cos(th(j))*psip(j)*L(j); 
nz(j)=V(j)/g*thp(j)*L(j)+cos(th(j)) ; 

alphag(j)=180/pi*sqrt(ny(j) A 2+nz(j) A 2)*m_i*g/(ro*V(j) A 2*Sref/2)/13; 
time(j)=time(j-1)+dtime; 

end 

tgo=time(end)-time(1); 

for i=l:Npt 

nX(i)=norm(X(1:3 , i)); 
nXp(i)=norm(Xp(1:3 , i) ) ; 
nXpp(i)=norm(Xpp(1:3,i)); 

end 

trys=[X' nX' Xp' nXp 1 Xpp' nXpp']; 

path=[time' X(1:3,:)’ V' th' psi' nx' ny' nz' tau' L']; 

%% Computing cost and penalties 
V_f=V(end); 

xdt=BC(14:16)+BC(17:19)*tgo; 
xdf=[V_f*cos(thf)^cos(psif); 

V_f*cos(thf)^sin(psif) ; 

V_f^sin(thf) ] ; 

nmax=40+(40-10)/(0-50000)*(X(3,:)-Re); 
if nmax<l A nmax=l; end 

J=[tgo; 

abs(dot(xdf , xdt))/norm(xdf)/norm(xdt)]; 

Py=max([0 , abs(ny)-nmax]); 

Pz=max([0 , abs(nz)-nmax]); 

Dtgo=tgoold-tgo; 

Alt fine=max([0,X(3,end)-Re-60000]) . A 2+min([0,X(3,end)-Re-9000]) . A 2; 
cost=norm([1,weight90]*J)+10*(Py A 2+Pz A 2)+100*Dtgo A 2+0.l*Altfine; 
costs(qq,:)=[cost;tau_f;J;Py;Pz;Dtgo;Altfine]; 

%% Animating iterations 


figure(100) , % set(gcf , ’Color ! , ’w’); 
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subplot(3,6,[1 2 7 8]) 

plot3(Pos_BR(1:300,1)/10 A 3,Pos_BR(1:300,2)/10 A 3,(Pos_BR(l:30 0,3)- 
Re)/10 A 3, ... 

' - . r ' , 'LineWidth' , 2) 

hold on, ylim([0 60]); %axis equal 

plot3(Pos_BR(14 0,1)/10 A 3,P o s_BR(140,2)/10 A 3, (Pos_BR(140,3)-Re)/10 A 3, . . . 

' A g' ,' LineWidth' , 2) 

plot3(Pos_SM(1:19,1)/10 A 3,Pos_SM(1:19,2)/10 A 3,(Pos_SM(1:19,3)- 
Re)/10 A 3, ... 

' . k ' , ' LineWidth' , 2) 

plot3(X(1,:)/10 A 3,X(2,:)/10 A 3,(X(3,:)-Re)/10 A 3, 'Linewidth' ,3); grid on 
plot3 (xf(1)/10 A 3,xf(2)/10 A 3, (xf(3)-Re)/10 A 3, f pk f , 'MarkerSize ',11) 
plot3(X(1,end)/10 A 3,X(2,end)/10 A 3, (X(3,end)- 
Re)/10 A 3, 'pr’, ’MarkerSize' , 9) ; 
view(-102,8) 

hl=legend( 'BM trajectory' , 'BM detection ',' Unguided ascend',... 

'Guided flight', 'Predicted intercept point', 'Actual intercept 
point' ,... 

'Location' , 'North' ); set(hi, 'FontSize' ,7) 
hold off 
%{ 

axis([0 1.Ie5 -le4 1.Ie5 0 7e4]/10 A 3) 

xlabel('Northing, x (km) ' ),ylabel('Westing, y (km) ' ),zlabel('Altitude 
(km)') 

subplot(3,6,[13 14]) 

plot (path(:,1),path(:,9), '—b', 'Linewidth',2) 
hold on; grid on 

plot (path(:,1),path(:,10), '-.g', 'Linewidth',2) 
plot(path(:,1),nmax(:),'r','Linewidth',2) 
plot(path(:,1),-nmax(:),'r','Linewidth',2) 
axis([10 time(Npt) -45 45]); 

xlabel('Time (s)'), ylabel('Load Factor (g)') 
leg2=legend('n_y','n_z','Dynamic 
constraints','Location','Best'); 

set(leg2,'FontSize',7) 
hold off 
subplot(3,6,3) 

plot(time,X(1,:)/10 A 3, 'Linewidth',2) ; grid on 
axis([10 time(Npt) Ie4/10 A 3 I.le5/10 A 3]) 
title('x_l (km)') 
subplot (3,6,9) 

plot(time,X(2,:)/10 A 3,'Linewidth',2); grid on 
axis([10 time(Npt) -Ie4/10 A 3 I.le5/10 A 3]) 
title('x_2 (km)') 
subplot(3,6,15) 

plot(time, (X(3, :)-Re)/10 A 3, 'Linewidth',2); grid on 
axis([10 time(Npt) 0 7e4/10 A 3]) 
title('x_3 (km)'), xlabel('Time (s)') 

subplot(3,6,4) 

plot(time, Xp(1, :)/10 A 3, 'Linewidth',2); grid on 
axis([10 time(Npt) -45 45]); axis 'auto y' 
title('x_l' ' /10 A 3 ' ) 
subplot (3,6,10) 
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plot(time,Xp(2, :)/10 A 3, 'Linewidth', 2); grid on 
axis([10 time(Npt) -45 45]); axis ’auto y' 
title('x_2’ ' /10 A 3 ' ) 
subplot(3,6,16) 

plot(time,Xp(3,:)/10 A 3,'Linewidth',2); grid on 
axis([10 time(Npt) -45 45]); axis 'auto y' 
title( 1 x_3' 1 /10 A 3') , xlabel('Time (s) f ) 
subplot (3,6,5) 

plot(time,Xpp(1, :)/10 A 2, 'Linewidth',2); grid on 
axis ([10 time(Npt) -45 45]); axis 'auto y' 

title('x_1' ' ' ' /10 A 2 ' ) 

subplot (3,6,11) 

plot(time,Xpp(2, :)/10 A 2,'Linewidth',2); grid on 
axis ([10 time(Npt) -45 45]); axis 'auto y' 
title('x_2' ' ' ' /10 A 2' ) 
subplot(3,6,17) 

plot(time,Xpp(3, :)/10 A 2, 'Linewidth',2); grid on 
axis([10 time(Npt) -45 45]); axis 'auto y' 
title('x_3 ' ' ' ' /10 A 2'), xlabel('Time (s)') 
subplot(3,6,6) 
hold on 

plot(qq,time(Npt)-time(1),'+c','Linewidth',1); grid on 
hold off 

axis([l 200 40 60]); axis 'auto y' 
title('Time-to-go (s)') 
subplot(3,6, [12 18] ) 
hold on 

plot (qq,cost,'+m','Linewidth', 1) ; grid on 
axis([l 200 0 100]); axis 'auto y' 
hold off 

title('Performance Index'), xlabel('Iteration') 

%} 

return 


D. COMMON FUNCTION - STANDARD ATMOSPHERE 
1. STatmos.m 


function [Density, Pressure, Temperature]^STatmos(alt) 

% Calculation of the 1976 standard atmosphere up to 86 km 
% Code source: http://www.pdas.com/atmos.htm 

% Run ezplot('STatmos',[0,86000]) to see the plot of density vs 
altitude 

o 

o 

% Author: Yakimenko, Oleg A. 

% Date: September, 27 2005 

% E-mail: oayakime@nps.edu 

o, 

alt=alt/1000; % Convert altitude from m to km 

%% - Initialize values for 1976 atmosphere 
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REARTH=6369.0; % Earth radius (km), depends on Latitude 

GMR= 34.163195; % Gas Constant 

htab=[0.0, 11.0, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852]; % Geometric alt 
ttab=[288.15, 216.65, 216.65, 228.65, 270.65,... % Temperature 

270.65, 214.65, 186.946]; 

ptab=[1.0, 2.233611E-1, 5.403295E-2, 8.5666784E-3, ... % Relative pres 

1.0945601E-3, 6.6063531E-4, 3.9046834E-5, 3.68501E-6]; 
gtab=[-6.5, 0.0, 1.0, 2.8, 0.0, -2.8, -2.0, 0.0]; % Temp gradient 

P0=101325.0; Ro0=1.225; 

%%- Convert geometric to geopotential altitude 

if alt>250 

alt = 100 

end 

h=alt*REARTH/(alt+REARTH) ; 

%% - Binary search for altitude interval 

i = i; 

j = 8; 

while j > i+1 

k=fix((i+j)/2); 
if h<htab(k); 

j = k; 

else 
i = k; 
end 

end 

%% - Calculate local temperature 

tgrad = gtab(i); 

tbase = ttab(i); 

deltah = h - htab(i); 

tlocal = tbase + tgrad*deltah; 

theta = tlocal/ttab(1); 

%% - Calculate local pressure 

if (tgrad == 0.0) 

delta = ptab(i)*exp(-GMR*deltah/tbase); % Isothermal layers 

else 

delta = ptab(i)*(tbase/tlocal) A (GMR/tgrad); % Non-isothermal 

layers 
end 

%% - Calculate local density 

sigma = delta/theta; 

%% - Current atmosphere parameters corresponding to Altitude alt 

Temperature=tlocal; Density=Ro0*sigma; Pressure=P0*delta; 
return 
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APPENDIX E: DETAILED DESCRIPTION OF TRAJECTORY 
SHAPING GUIDANCE (REPRODUCED FROM [10]) 


A method that overcomes fatal characteristics of modem guidance laws is the key 
driving factor in the development of this advanced guidance law. The guidance law will 
determine the near-optimal flight path from the interceptor position to a predicted target 
position for the interceptor to follow to intercept, and then derive the set of control 
commands necessary to execute that flight. This section will describe a method for 
deriving that trajectory using calculus of variations based on three cues: high-order 
polynomials as a reference function for the flight path, a preset thrust history as one of 
the controls, and a few optimization parameters. The trajectory optimization problem is 
then converted into a nonlinear programming problem and solved numerically. 

A. PROBLEM STATEMENT 


Among all admissible trajectories, 


that satisfy: 

z(t) = {z l (t),z 2 (t),...,z r (t)} T eS 

S = jf(0e2 r c£ r }, te[t 0 ,t f ] 

(4.A.1) 

1. 

The system of differential equations (dynamic constraints): 



Zi = f i (t,z,u,a), i = 1 ,r 

where the vector of controls is 

(4.A.2) 


u(t) = {u l (t),u 2 (t),...,u m (t)Y, m<r, ueU m (^E m and the 

missile parameters is a = {a v a 2 ,...,a p ), a e. A p c= E p ; 

vector of 

2. 

The initial conditions: 



z(t 0 )eS 0 , ^=jz 0 eZ r cT) 

(4.A.3) 


u(t 0 )eR,, R 0 ={u 0 eU m ^E m ] 
and the final conditions: 

(4.A.4) 


z(//)e5/, S, = (z,EZ'cr; G 1 (z(l l )) = 0,j = \J^ 

(4.A.5) 


m 

II 

m 

a 

n 

fra 

3 

(4.A.6) 

3. 

The constraints imposed on the state space 
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(4.A.7) 

on the controls 


I (t, Z, u) = {4j (t, Z, u), 4 (t, Z, U), ..., 4 (t, Z, u)Y > 0 

(4.A.8) 

and on the controls derivatives 


ft{t,u,u) = |^j(t ,U ,u), 7T 2 (t ,U ,u), ..., 71^(t ,U ,l/)| > 0 

(4.A.9) 


Find the optimal trajectory, z t (t ), that minimizes the integral function 

‘r 

J = K(x 0 ,x f ) +1 L(t,z,u)dt (4.A.10) 

*0 

and the corresponding optimal controls, u ,{t ), where K, L are defined functions. 

B. CALCULUS OF VARIATIONS 

Calculus of variations deals with functions of functions, termed functional, 
instead of functions of some variable or variables as in ordinary calculus. Specific 
interest is in the externals of these functionals - those making the functional attain a 
maximum or minimum value [Ref 15]. There are two broad categorizations of methods 
to solve these problems, indirect methods and direct methods. 

Indirect methods resolve the problem into a differential equation, usually via the 
difference of a series of equations of motion and thus solve the general theory of partial 
differential equations. This method does not assume anything about the solution, leading 
to a very complete and perfectly precise answer; however, one must integrate the 
resultant equation to derive the extremals. The precision greatly increases the 
computational complexity, and therefore the time required to solve, for negligible gain in 
optimality. Further, these differential equations are difficult to integrate except in the 
simplest of cases, and nearly impossible to program [Ref 19]. This approach is further 
complicated by the need to solve the problem in a specified fixed region instead of in the 
small neighborhood of some point. These difficulties can be overcome by using direct 
methods, which do not reduce the variational problems to ones involving differential 
equations. 
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The fundamental idea of direct methods is to consider a variational problem as a 
limit problem of the extreme of a function of a finite number of variables that can be 
solved by numerical methods. Two basic direct methods are the Rayleigh-Ritz and 
Galerkin methods, which assume the solution to be an unknown function but of a certain 
form containing a set of unknown coefficients (themselves functions of the boundary 
conditions), which are then found by minimization. The practical result is that the 
problem has been reduced from calculus to algebra, but at the cost of a significant 
increase in the number of simultaneous equations to be solved. It is for this reason that 
little work was done in this field until the advent of computer technology. The possible 
solution set is restricted to a smaller space than the original equation because of the initial 
assumption regarding the solution, but the resultant problem can be programmed and 
quickly solved using computers; however, the solutions are only an approximation of the 
original solution. Therefore these solutions can only be regarded as near-optimal 
solutions. 

Professor Taranenko [Ref 19] first applied the ideas of direct methods and the 
combination of Ritz and Galerkin methods to the problems of flight dynamics by 
identifying a reference function for the flight vehicle’s motion and velocity 

T = x w + ( x i f ~ */o ) ~ r ° + < h,(r), / = 1,4 (2.1) 

T f ~ T o 

where x,, x 2 , x,, are the Cartesian coordinates for the flight path, x 4 is the velocity, and 
0 ( (r) is a continuous, unequivocal, differentiable function satisfying the boundary 
conditions O ( (r 0 ) = 0.(r f ) = 0. Any function that satisfies those conditions would be 

acceptable for use. Taranenko called r a virtual arc. It is this critical variable that allows 
the separation of the spatial trajectory from the velocity and thus optimize one or the 
other, or both independently. The specific task determines the appropriate choice of r, 
but in general any continuous, monatomic function is acceptable: time, path, energy, etc. 
The remaining state parameters and flight controls are then determined by solving the 
inverse problem of flight dynamics. Rather than starting with the control time histories 
and integrating them to determine the flight path, this method starts with a flight path and 
determines the control time histories necessary to create it. 
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In order to implement a solution to this problem that can be solved in real time, a 
further restriction must be made. The continuous problem must be discretized to reduce 
the infinite variational problem to one of optimization of few parameters at numerous 
sampling points. This allows for the optimization of the planar trajectory of the missile at 
several points along the path by presetting the state variables and one of the controls’ 
time histories, and then solving the inverse flight dynamics problem. 

These methods have all previously been used for off-line optimization of 
trajectories, but none has yet been applied to the real-time onboard optimization of a 
missile flight path. Yakimenko detailed a method called a Direct Method for Rapid 
Prototyping, which he applied to short term spatial trajectories of aircraft maneuvers 
using fixed boundary points [Ref 19]. The developed program presented here uses similar 
numerical method to provide a near-optimal spatial trajectory that is completely defined 
by a few optimization parameters, but as will be discussed shortly, has fluid final 
boundary conditions. 

Though the method artificially limits the possible trajectory variations, it does 
guarantee the following [Ref 19]: 

1. The boundary conditions are satisfied a priori, 

2. The control commands are physically realizable and smooth, 

3. Only a few variable parameters are used, thus ensuring that the iterative 
process converges well, 

4. The near-optimal solution is very close to the optimal one. 


A simple two dimensional variation program of a similar 7 th order system will 
demonstrate how the direct method varies the flight path according to the boundary 
conditions. The boundary conditions are 


o 

II 

o 

J-T 

o 

II 

o 

X lf Z 

= 1 

x 2f =l 

Xj ' 0 = 0.2 

x — ] 

A 20 1 

x[f = 

0.1 

x' 2f = -1 

4 = 0.1 

II 

O 

x” = 

0.1 

x” 2f = 0.1 

x"\ 0 = var 

x'" 20 = 0.1 

x m = 

1/ 

0.1 

x"' 2f = 0.1 
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where the third derivative of the initial Xj condition has been set to vary according to 

x"\ 0 = {-0.4;-0.1; 0.2; 0.5} 
and the length of the virtual arc varies according to 

T f = 1,2,...,10 

The resulting set of paths is shown in figure 26 and the first derivative of the path is 
shown in figure 27. It is important to note that the first derivative is not “velocity” but is 
instead the “rate of change of the path”. It is proportional but not equal to velocity, 
specifically because of the virtual variable z as discussed previously. Each line 
represents a different choice of z f , showing that by varying that value the length of the 

path can change drastically. The algorithm developed here will use this technique to 
derive the optimal flight path. 


d 3 x/dx 3 | 0 =-0.4 



c\j 


0 0.5 1 1.5 2 


d 3 x/d T 3 | 0 =0.2 
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d 3 x/dx 3 | Q =0.5 



d 3 x/dx 3 | 0 =-0.1 



Figure 37. Variation of Path with z f and x"' (after [Ref 18]) 
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Figure 38. Variation of First Derivative of Path with r f and x'" 0 


C. PROGRAM DEVELOPMENT 

1. Boundary Conditions 

Using equations (3.2) and (3.3), in addition to the controls n x , n y ,n _, it is possible 
to construct the vector of state variables z = {x l ,x 2 ,x 3 ,V,0,¥} T and the vector of controls 
u = {n x ,n ,n z } T . The model for thrust, drag, and missile characteristics is the same as 
previously used in Chapter 3. 

The beginning assumption is that the following data is known by the interceptor’s 
onboard computer: 
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Interceptor 

Target 

Body Frame 

X l (t 0 ) = x l0 

^2 (^0 ) — ^20 
(^0 ) — A$0 _ 

V(t 0 ) = V 0 -^10 

m) = ^0 ^ *20 *10 

&( t o) = &0 *30 ^*20 

n y(t*) = n y« *30 

n z (t 0 ) = n z0 J 

Xi(t / ) = x 1/ 

x 2 (t f ) = x 2f 

X 3 (t f ) = X 3/ 

'P(t / ) = 'P / ^ 2/ 

9(tf) = df J x 3/ 

Earth Centered 
Inertial 

xT (o, vf (t) 

y M (t), if (t), ifft) 
xf(t), xf (0, xf (t) 

xf (0, xf (0, xf(0 
xf(t), if (0, if (0 


Table 9. Interceptor Known Data 


The target data will be used to determine the final boundary conditions of the 
interceptor missile according to the time until intercept, At: 

Position: 


SM 


BR 


X \ f = X 10 + X 10 ^ 


SM 


BR 


^ 2 ^ ^20 I JC20 / 


SM 


BR 


x 3 f — x 20 + x 20 At 


Heading and Flight Path Angle: 


qSM QBR ^ w ^ ere QB& _ arct g 


71 


x 3 r 


Ri 


BR2 •BR2 
I *^2 


•BR 


X, 


(//f = yr B f R , where y/ B f R =arctg- 


x, 


(4.C.1) 


(4.C.2) 


which, when combined with reasonable estimates of the final values of the velocity, V, 
and the time rate of change of velocity, V , heading ('P = 0), and flight path angle 
(<9 = 0), yields the final conditions: 
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(4.C.3) 


and 


xf? = V f cos 6j cos y/ f 
x“ = V f cos 0 f sin y/ f 
x 3 ™ = V f s i n 0 f 


xff = V f cos Q f cos f - OjV f sin 0 f cos *¥ f - V F / F / cos 6 f sin *¥ f 

xfj = V f cos 6 f sin W f - 6 f V f sin 6 f sin m f + 'i' f V f cos 6 f cos V F / (4.C.4) 

* 3 / = ^ / s i n + QfVf cos 0 f 

In order to ensure a smooth flight path at the initial and final points, an additional 
constraint of 

■■■SM _ ■■■SM _ n 

-h; — A i / — u 

x™ = X™ = 0 (4.C.5) 

■■■SM _ ■■■SM _ n 

A 3 

will be imposed on the system at the initial and final conditions. 


2. Separating and Recombining Space and Time 


As discussed previously, in order to independently optimize the spatial trajectory 
and the velocity, the reference function will be derived as a function of z . The boundary 
conditions cannot, therefore, be defined as functions of time derivatives as in equations 
(4.C.3) - (4.C.5). A connection between the spatial and time domains must therefore be 
introduced, X , which is defined as 



(4.C.6) 


and is termed the virtual speed [Ref 19]. This allows for the independent variation of the 
speed profile along the same paths according to any other convenient reference. In this 
case the known thrust profile,^, will be utilized by integrating the third equation of 
(2.B.3) and applying the virtual speed 


t/v \ / • n\ dT g(n x -smff) 

V ( r ) = g( n x ~ sin 0) — = -- 

at X(z) 


(4.C.7) 
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Further use of the virtual speed allows the recalculation of the initial and final boundary 
conditions, transforming them from the time frame to the spatial frame. For this, the 
obvious relations 


X i (T)=^ L ^ = x' i (T)A(T) 
dr dt 


x,(r) = 


d(xi(r)A(T)) dr „, 2 , . 


(4.C.8) 


= x"A + x t A' i = 1,2,3 


dr dt 

which when rearranged defines the first and second derivatives of the missile coordinates 
as 

x,' = r X,. x';= X- 1 [x. -x t A'] I =1,2,3 (4.C.9) 

Using the values of A and A' defined as 


A = V A' = V V A = V A' = V V 

/L 0 V 0 5 /L, o y o y 0 /l f y f’ /l f y f y 1 


-1 


(4.C.10) 


3. Reference Trajectory 


The knowledge of the initial and final position plus the initial and final conditions 
of the first and second time derivates allows for the construction of a 7 th -order 
polynomial (the maximum orders of the time derivatives of the missile coordinates at the 
initial and final points plus one) to describe the reference function of the aircraft 
coordinates x ; (i = 1,2,3). The following are introduced as the reference functions [Ref 


19]: 


Or written another way, 


x 




k =0 
5 


(max(l,&-2))!r* 

~k\ 


(max(l,A:-2))!r 
M = -(tTl)! 


k=l 

5 


x%r) = Y J a ik r 


k -2 


k =2 


x ( W = Z(*-2Kr‘- J 

k=3 


(4.C.11) 
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— G^ + + G^T + G^T^ + G^T 

x "i (?) = a i2 + a n z +1 a j4 r 2 +1 a i5 r 3 + ^ a i6 r 4 +1 a n z 5 

x'.(r) = a a +a i2 r + -a i3 r 2 +-a i4 r 3 + — a i5 r 4 + — a i6 r 5 + — a n T 6 
2 6 12 20 30 

, , 1 2 1 3 1 4 1 5 1 

x,(T) = a i0 +a a T + -a i2 T +~a i3 r + — a iA r + — a i5 r + — 


(4.C.12) 


a i6 r + 


210 


a n T 


The coefficients can be determined by solving the equations simultaneously, 

V.. 3 


n 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 
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0 
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0 

0 

0 
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1 

0 

0 

0 

0 
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z f 
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2 L f 

It 3 

6 V 

J_7- 4 

24 if 

ir 5 

60 if 

J_7- 6 

120 v 

210 r / 

0 

1 

T f 

It 2 

2 7 

l r 3 

6 7 

if 

12 V 

20 if 

if 

30 if 

0 

0 

1 

r / 

l r 2 

2 V 

if 

3 if 

if 

47 

if 

1° 

0 

0 

1 

r / 

7 

V 

f 


a to 

a,. 


a i 2 


°i 4 


7 6 


f 


f x 3 

X U) 


*/o 

ft 

X i0 


X 


Xu 


x , 

x \ 

x" 


(4.C.13) 


if J 


Finally, by substituting the corresponding values of x i0 , x ' 0 , x” {) (i = 1,2,3) for 
r 0 =0, and x if ,x\ f ,x" (7 = 1,2,3)for r = r f (where r f is the first optimization parameter, 

the virtual arc), results in a set of 24 linear algebraic equations for 21 unknown 
coefficients a ik (i = 1,2,3, A: = 0,1,..,7) 


U i0 X i0 


a i\ = X i0 


a j2 x i0 


a n = x 


-Ax’" +16x,"; 60x,;-120x" 0 -360x' -480x,' 0 840x, r -840x i0 

a i4= - Z -+-2-+- H -+- ~A - 


'/ 


V 


'/ 


3Ox!! + 60x;;; -420x" + 600x,; 2340x' - 2700x; o 

a i5 = - ~ -+- — -+- —A -+ ■ 


V 


v 


v 


v 

-5040x (/ +5040x i0 


v 


(4.C.14) 


-60x! - 80x™ 780x! -900x" 0 -4080x7 -4320x' 0 8400x, r - 8400x, 0 

// — ---1---1---1--- 

U i 6 _3 ' _4 ' _5 ^ _6 


'/ 


V 


'/ 


35x" + 35x"' -420x" + 420x" 0 21004+2100x1 -4200x, + 4200x f0 

a n =--4- + -b- + --6- + -b- 


V 


V 


V 


V 
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4. Inverse Dynamics 

The numerical solution develops a reference trajectory over a fixed set of N points 
equidistantly placed along the virtual arc. The virtual interval is 


At = —— 
N -1 


(4.C.15) 


and the corresponding time interval is 


f 3 


A tj = 2 


2X/1> 

V 1 _ J 

V J +V J -1 




(4.C.16) 


where V. and V J _ ] is determined by integrating equation (4.C.7). N is any convenient 
number, chosen to be 100 in this program. 

The value of V, also allows for the determination of both 6 and T* by rearranging 
equations (2.B.2) and substituting the virtual velocity 

r 

X. . 

(9 = arctg- 

Jxl.S+xi/ 

(4.C.17) 




x. 


¥j= arctg — 


2 ;j 


X Uj 


The controls n y and n_ are found by rearranging equations (2.B.3) to be 

V. . 

riy.j = — x ¥ j cos 6j 
§ 

n =—(&j+ cos6») 

g 

where the angular derivatives are determined as 

Q _ cos 2 Q X 3 "( X l ' 2 + X 2 2 ) ~ X 3 ( X l"+ X 2 ) 

1 (xi 2 +x' 2 2 r 2 

- U/ X I X 2 ~ X l X 2 ) 


T = cos T 1: 


(4.C.18) 


(4.C.19) 


x, 
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5. Cost and Penalty Functions 

Finally, the calculation of the flight path results in a set of functions that must be 
minimized, which occurs through a Cost Function (CF) and a Penalty Function (PF). 
These functions must be carefully chosen to ensure that the optimal path is truly feasible, 
desirable, and obtainable. A simple example of a CF is J = t f , the minimal time 

problem, or J = \u(t )\, the minimum fuel problem, though there is no limit to the number 
or variation of the Cost Function. 


In this case, the CF was chosen to optimize three properties simultaneously; 
minimize the length of the virtual arc, r f , minimize the time to intercept, t , and 


maximize the impact angle of the interception. The CF for this program is written as 


J = w t k t r f + w 2 k 2 1 + w 3 k 3 


V BR .y 


SM 




(4.C.20) 


Each item must be scaled appropriately using the scaling factors k v k 2 ,k 3 , so that 


they are all roughly equivalent when optimized, e.g. the anticipated intercept time is 
counted in tens of seconds while the cosine of the impact angle will vary from zero to 
one. Failure to weight them properly will skew the results of the cost function. 
Additionally, through the weighting functions w t ,w 2 ,w 3 , a trade-off analysis can be 


conducted and variables can be included or excluded as desired. 


The first two variables, r and t , are necessary to ensure the system’s optimal 

solution is actually physically realizable. Without them, the system will continue to 
optimize the intercept well beyond the capabilities of the missile or even physical reality. 
One example is that, in a purely mathematical sense, there is no problem with a negative 
velocity (yet in the physical world that makes no sense) and the program might try a 
program that would intercept after the missile velocity has gone past zero and into 
negative numbers (of course, the missile would stop flying long before it even reaches 
zero). One might be tempted to include a myriad of parameters to cover all possible 
eventualities; however, including these two parameters sufficiently accounts for nearly 
every physical limitation and eliminates the need for having a long list of cost variables. 
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The third variable in the CF was chosen in order to maximize the angle of impact 
upon interception. This reflects the need, described in Chapter 1, to maximize the kinetic 
energy in order to ensure the interceptor disables the target. The cost is calculated using 
a simple dot product relationship 


cos 0 = 


7 i/-* SM ""•* BR \ 

dot(x f , x f ) 



ZSM 

Z BR 



Xf 

x f 



(4.C.21) 


which will have a minimum value of zero when the impact angle is maximum. 


The PF is chosen to ensure that certain conditions are not violated or exceeded, 
such as physical limitations. In this case, the PF is on the maximum acceleration in the y 
and z direction. Zarchan showed that acceleration capability is dependent on altitude and 
speed [Ref 21], The PF has been set up to reflect this by varying between 40 g’s at sea 
level to 10 g’s at 50,000 ft. 


These are not the only CF or PF variables that can be included. The choice of 
variables to include is situation dependant and may be modified to meet whatever the 
needs of the situation demand. 
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